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Deflection and Stress Patterns in Plate-Type 
Structures, Theory vs. Test Results 


NATHANIEL J. HESS* 
General Electric Company 


Summary 


This report presents a practical analysis of the deflection and 
stress distribution throughout a plate-type structure with vari- 
able stiffness, loading, and boundary conditions. Its application 
is to aircraft and missile wing and fin structural design, and the 
determination of the deflected shape of an airfoil surface under 
aerodynamic loading. The basic assumption used for the ap- 
proximation of the rigorous theory is that of idealizing deflec- 
tions. In this manner, the chordwise deflection at any spanwise 
station is expressed in the form of a parabolic equation, where the 
coefficients are unknown functions of the spanwise stations. By 
utilizing the idealization, each spanwise station has three degrees 
of freedom that must be determined. This analysis allows the 
theoretical matrix plate equations to be programed easily on a 
high-speed digital computer; once programed, a practical wing 
may be designed by two or three stress analysts in a short period 
of time. 

The particular theory developed in this paper is for application 
to full-depth honeycomb construction, which is very efficient 
weightwise for a thin low-aspect-ratio airfoil surface. The 
problems presented are for airfoil surfaces attached to the side of 
a shell-type structure with full and partial root support. Com- 
parisons of root support are shown in order to demonstrate the 
variation of deflection and stress distribution. 

In order to establish the validity of this analysis, and to verify 
that the results obtained may be applied to a plate-type struc- 
ture, a test program was undertaken. The theory was then ap- 
plied to two typical problems, and an experimental program was 


conducted. The results were then compared. 
Symbols 
{] = square or rectangular matrix 
{} = column matrix 
[|]? = matrix transpose 
[]~ = matrix inverse 
[A] = stiffness matrix 
la] = stiffness submatrix 
[B] = geometry matrix 
¢ = length of root chord 
[D}] = differentiating matrix 
D(x, y) = Plate stiffness Et?/[12(1 — u?)] 
E = modulus of elasticity 
é = station spacing (constant) 
$ = (2)th station 
L = span length 
| M,,} = moment distribution about the x axis 
N = number of stations 
{P} = loading matrix 
t = plate thickness (total) 
U = strain energy 





Paper presented at the IAS National Summer Meeting, Los 
Angeles, June 28, 1960. Received February 25, 1960. Revised 
and received February 24, 1961. 

* Functional Engineer in Space Structures Engineering, Missile 
and Space Vehicle Department. Formerly Principal Missiles 
Systems Engineer in charge of design and analysis of aerodynamic 
surfaces at Republic Aviation Corporation. 
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V = potential energy 

fw} deflection matrix 

x = variable distance along the chord 

y = variable distance along the span 

{33 = functions of deflection (7 = 0, 1, and 2) 
u = Poisson’s ratio (constant) 

o = principal stress 

Or = chordwise stress 

Cy = spanwise stress 

Try = shear stress 


Introduction 


AS AIRCRAFT FLIGHT SPEEDS increase into the super- 
sonic and hypersonic range, and weight becomes 
of increasing importance, it has become necessary to re- 
fine the structural analysis. Due to the increased 
speed of aircraft and to minimize drag, the aerody- 
namic surfaces have become thinner, the aspect ratio 
has become lower, and the leading edges have become 
swept back. The conventional beam theory, as applied 
to high-aspect-ratio surfaces, is no longer adequate, and 
its use could result in an error in the stress distribution 
of a few hundred percent. In beam theory, bending 
and torsion are uncoupled due to the existence of an 
elastic axis; however, in the problems presented in the 
following discussion, this is not the case. Shear and 
moment diagrams are meaningless in plate-type struc- 
tures as beam deflections are not simulated. This makes 
it necessary to employ a more accurate analysis when 
approaching a low-aspect-ratio surface or that of a delta 
wing which is described in this paper. 

This analysis will be useful in establishing new tech- 
niques for future problems, where similar plate-type 
structures are encountered. Examples of typical prob- 
lems are as follows: (1) triangular plates simply sup- 
ported or fixed on two sides; (2) building design where 
large roof overhangs exist; and (3) long, wide bridge 
spans, etc. In these problems, where the idealization 
of deflections represents the accuracy of the analysis, 
algebraic equations may not be representative of the 
vertical displacements. In some cases, a trigonometric 
series may be more desirable. 

The importance of this paper is in its presentation of 
practical solutions to theoretical problems. In this 
manner, the nonspecialist may readily obtain a “‘sense 
of feel” for the deflection and stress distribution for 
these and similar problems. Once the distribution is 
understood, a more efficient structure can be designed. 

For the purpose of comparing theoretical and test re- 
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sults, a delta planform plate with a 60° swept leading 
edge was used. The geometry was based on a particular 
case which was to be used on a practical wing design at 
a later date. However, the analysis method is not 
limited to a delta planform; it may be applied to any 
plate meeting the basic boundary conditions and de- 
flected-shape characteristics. 

Various parameters were imposed on the wing con- 
figuration so that a simple experimental model could be 
constructed; they were as follows. 

(1) A 60° delta plate was used first with full, and 


then partial, root support. 
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(2) A uniform thickness was used throughout. 

(3) General unit loadings were applied. (In this 
paper, only unit loadings at the apex are shown; how- 
ever, other unit loadings at arbitrary points were 
checked with similar results.) 

(4) A uniform load of 1 psi was applied to the fully- 
fixed plate by means of the theoretical analysis, and a 
comparison was made with NACA theory and test 


results. ! 


Method of Analysis 


The structural analysis presented is for an elastic, 
homogeneous, and isotropic plate, where deflections 
may be large ccmpared to the plate thickness. The 
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total potential, which is the difference between the 
strain energy in the plate and the potential energy of 
the external loads, is minimized with respect to the 
unknown redundants of the idealized deflection. 

The overall differential equation, as a function of the 
deflection for the bending and twisting of a plate, repre- 
sents the strain energy and is given by the following 


Wy? 
) + 


equation: 


9 be 1 ct x) me ' (<= Pe 
~— 2J0 J hix) me Ft ‘ Ox? 


Oy 
2(1 | (S5.) ae | ixdy (1) 
2(1 — u _ > dxdy 
Oxdy Ox? Oy? ee 
The energy of the external loads is given by 
L f(x) 
V=- f f P(x, y) W(x, y)dxdy (2) 
0 Si(x) 


Eqs. (1) and (2) are derived by Timoshenko.’ 

The method proposed for the solution of these equa- 
tions is to assume a deflected shape for the chordwise 
stations. The equation used will be that of a parabolic 
shape as recommended by Stein, Anderson, and Hedge- 
peth! and is given by 


W = do(y) + xor(y) + x°¢o(y) (3) 


The coordinate axis and their direction are shown in 
Fig. 1. For each chordwise station (a function of y) an 
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Theoretical deflection pattern of a partially fixed plate 
with a 70-lb load applied at the apex. 
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initial deflection at the y axis is given by @po, the slope 
of the twist by ¢:, and the chordwise bending by ¢». 

The basic derivation for the deflection matrix equa- 
tion is presented by Stein and Sanders in reference 3. 
In this reference, deflections are determined theoreti- 
cally for a thin, low-aspect-ratio wing of built-up con- 
struction. Two problems are presented, one being a 
plate with a symmetrical clamp and loading, and the 
other having simple supports and antisymmetric load- 
ing. The problems presented in this paper are for plates 
with full and partial fixity along one edge with asym- 
metrical loading. Although the basic boundary con- 
ditions in reference 3 had to be modified for this prob- 
lem by changing the stiffness and geometric matrices, 
they have been invaluable for their theoretical matrix 
development used in this paper. 

The idealized deflection represented by Eq. (3) is then 


THE EQUATIONS OF THE ANALYSIS 
ARE NOT APPLICABLE IN THIS 
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substituted into the energy expressions denoted by 
Eqs. (1) and (2). All integrations are done numerically 
by the trapezoidal rule, and differentiating is done by 
finite differences. 

The expression is then minimized with respect to the 
redundant deflection parameters to obtain the following 
relationship: 

a(U + V)/d#i, = 0 (4) 
and OU /0¢i; = —OV/0¢;; (5) 
This expression yields 
[Alldu} = [Bt P} (6) 
lou} = [4] 1B] P} (7) 


where the redundant deflection parameters ¢ are given 
as a function of the stiffness, geometry, and the loads; 
but 


and 


{W} = [B]* {4x} (8) 
and, therefore, 
{W} = [B]"[A]-[B]{ P} (9) 


which is the basic deflection equation used in this 
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analysis. In this manner, the geometric properties 
denoted by [8], stiffness parameters denoted by [A], 
and the loads denoted by | P}, are submitted to a pre- 
programed high-speed digital computer, and the deflec- 
tions are obtained directly. The rigorous development 
of Eq. (9) is presented in the Appendix to this paper. 

In order to determine the stress distributions, the 
following equations” are used: 

















Et = + ow (10) 
yy, = —— SS 
“ — 2(1 — u?) \ dy’ Ox? 
Et (or + od | (11) 
= u— 
7 2(1 — u?) \ dx? oy" 
Et(L — u) 2&W (19) 
fT, , = ze _ — 
vw 6201 — u? 2) dvdy 
o = (a; + 0,)/2 + 1/2V (0, — oy)? + 4try? (13) 
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Fic. 6. Deflection distribution along chord stations for a 1-psi 


uniform load on a fully fixed plate. 
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Once again Eq. (3) is substituted into Eqs. (10), (11), 
and (12), and after simplification, the spanwise stress, 
o, for example, reduces to: 
mes: |) = aay UDalleal + x1Dalton + 
x?[Ds]\ d2} + 2ue\¢(MOD)}] (14) 
Eq. (14) is developed in the Appendix. This equation 
is the basic expression for the spanwise stress dis- 
tribution throughout the plate. By knowing the deflec- 
— tion parameters from Eq. (9) and the plate stiffness 
D(x, y), the stresses at any point in the plate may be 
i wees determined. 
This analysis was then applied to the plates shown 
in Fig. 1 for various arbitrary loadings. The only dif- Pec. &. Test plete ba fir accents 
ference in the theory between the two plates was the 
method of determining the plate stiffness parameters. 
yperties These parameters are determined by integrating along 
ry [A], 2 and 3; theoretical stress patterns are shown in Figs. 
r | pre- ! and 5. 
deflec- A uniform load of | psi was applied to a fully fixed 
1200 zm Y ‘ a - Sn ‘ ‘ 
pment = a $ plate in the second problem. . This was done so as to 
paper. 7 check this analysis against Stein, Anderson, and Hedge- 
ns, the 3° peth.!. This analysis is the same as the first with the 
— ” exception of a modification to the load matrix. The 
= © ol oe deflection and stress distribution along various chord 
(10) => a FULL ROOT FIXITY stations are plotted for this problem and are shown in 
> 800} — THEORY Figs. 6 and 7. These are compared with the NACA 
i --- NACA THEORY test program and theoretical analysis. 
(11) © c © TEST (REF. 4) In order to compare the theory with the experimental 
oF program, Figs. 2, 3, 4, and 5 are replotted along the 
(12) 2+ chord stations, superimposing the test data. This is 
cs done in Figs. 12 through 17. The significance of these 
** figures is expressed under ‘‘Discussion of Results.” 
(13) 2 
Experimental Program 
An experimental test program was conducted so that 
Y=60% the stress and deflection pattern obtained from the 
o | —— | : lle theory could be verified and the difference between 
0 20 40 60 80 100 theoretical and actual values could be obtained. Pho- 
x/C (PERCENT OF ROOT CHORD) tographs of this test assembly are shown in Figs. 8 and 
Fic. 7. Stress distribution along chord stations for a 1-psi 
x uniform load on a fully fixed plate. 
the chord from one edge to the other. In the partially ir’ le ee 
fixed plate, the root attachment boundary was inte- : ; 
grated from 8.34 < x < 25.76, and in the fully fixed 
plate it was integrated from 0 < x <41.56. The neces- 
sary input parameters for Eqs. (9) and (14), as outlined 
in the Appendix, are calculated and submitted to digital 
programing. In the problems checked for this paper, 
it was found that o, and 7,, produced negligible stresses. 
For all practical purposes, the principal stress o was the 
same as o,. The stress shown in all the figures will 
i00 therefore be cy. 
The first load condition considers a 70-lb load ap- “ 
psi plied to the apex of the plates shown in Fig. 1. Theoreti- : = BS E: : 
cal deflection patterns for these plates are shown in Figs. Fic. 9. Strain recording equipment. 
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9. Fig. 8 displays the test plate in its jig fixture. The 
plate is simply supported along its centerline and sym- 
metrically loaded so as to insure that the deflection and 
spanwise slope along this line remain zero during load- 
ings. On one side of the plate are the deflection gages. 
The gages with the large faces have a travel of four- 
tenths of an inch and are graduated to one ten-thou- 
sandth of an inch. Those with the small faces have a 
travel of 1 in., except the gages at the apex which have 
a travel of 2 in., and are graduated to one-thousandth of 
an inch. 

On the other side of the plate are the strain gages. 
These gages are rosettes, and were placed parallel, per- 
pendicular, and at a 45° angle to the support. The 
gages parallel and at a 45° angle produced negligible 
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rer, I/4 IN. 2024-T4 
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FIXED SUPPORT 
48 IN 255 => ¢ 





41.568 IN 


B) PARTIAL ROOT FIXITY 
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i - 1/4 1N. 2024-T4 
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PINNED CONNECTION 
48 IN § —;—t¢ 
| “1 
‘3 les 5.01N 


8.34 IN—<—e——__»-|7.42 IN. 
+— 41.568 IN——» 


Fic. 10. Geometry of the test plates with full root and partial 
root fixity. Note: The pinned connections restrain the plate 
against vertical deflections, but do not resist bending moment 


stress readings. This indicated that only the gages per- 
pendicular to the root, indicating normal stresses, were 
important and that the principal axis was parallel to 
the support for all practical purposes. For the tests 
described in this report, the loads were placed at the 
apex of the plate symmetrically about the support. 
These loads were introduced into the plate by hanging 
calibrated 10-lb lead weights, as shown in Fig. 8. Fig. 9 
displays the strain-recording equipment which con- 
sists of two SR-4 strain indicators and two switching 
and balancing units. The figures presented in this 
paper were for loads of 70 lbs applied at the apex of the 
plate. To obtain this increment, a 10-lb load was 
placed at the tip, and zero readings were taken on the 
dial and strain gages. An additional 70 lbs were then 
applied at the tip and the data were reduced for the 
increment. 
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Location of the deflection and strain gages on the test 


plates 


Two problems have been selected in order to show 
the effect of the boundary support and how radically it 
can displace the stress and deflection pattern. Fig. 10 
shows the geometry of these models; Fig. 10a shows the 
plate fully supported to resist the shear and moment 
along its centerline, and Fig. 10b shows the plate cut 
back along the centerline so that only a partial root 
support is felt by the plate. In this way, the entire 
moment is transmitted through this portion of the plate, 
and the shear is carried by this portion and by the shear 
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for a 70-lb load applied at the apex of a partially fixed plate. 
Note: The theoretical curve was obtained by fairing in smooth 
spanwise curves between the chordwise distributions. 
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pins added at the ends of this support. Fig. 11 shows a 
plan view sketch of the test plates giving the respective 
locations of the deflection and strain gages. 

Because of interference with the jig support, strain 
gages could not be located along the centerline of the 
plate, and the stress concentrations inboard of the 4.2- 
percent chord station could not be measured. In the 
partially fixed plate, shown in Fig. 5, the stress con- 
centrations in the aft region along the centerline were of 
particular importance. To check out this area, a 
photostress technique was employed to obtain the 
photograph shown in Fig. 18. The photostress tech- 
nique, known as the Zandman method, is a quantitative 
stress analysis combining the principles of photoelas- 
ticity, strain gages, and brittle lacquers. The structure 
is coated with a double refraction or birefringent ma- 
terial in a critically stressed area which is then subjected 
to loading. The stressed area is then illuminated with 
polarized light, and a complete color picture of the 
stress distribution over the entire area can be observed 
through an optical viewer. The points of stress con- 
centration are highlighted by different colors, each 
color representing a different stress level. When 
photographing this pattern with a black and white 
camera, fringes appear at various stress levels, as shown 
in Fig. 18. 

The results obtained from the experimental models are 
plotted in Figs. 12 through 17. They are superimposed 
upon the results obtained from the theoretical analysis 
and are labeled ‘“‘test results.”’ 


Discussion of Results 


Comparisons of the theoretical analysis and the test 
results for the plates described in Fig. 1 are shown in 
Figs. 6, 7, and 12 through 17. Figs. 6 and 7 show the 
comparison of deflection and stress distribution on a 
fully fixed plate. The agreement shown in this problem 
gives the necessary condition that the analysis yield 
practical results of sufficient accuracy to render it 
adequate for design purposes. In practical wing design, 
instead of placing a uniform load on the surface, a load 
matrix representing an airload distribution is used. 
Figs. 12 and 13 show comparable results of theory and 
test for a tip load on a fully fixed plate. Again, theoreti- 
cal versus test results indicate that concentrated loads 
yield an accurate analysis on a fully fixed plate. On 
practical structures, unit analysis must be calculated 
for concentrated loads such as pylon and elevon reac- 
Fig. 14 displays a comparison of deflections of a 
Figs. 15, 16, and 17 


tions. 
tip-loaded partially fixed plate. 
show the comparison of the theoretical results and test 
results for the stress distribution on the partially fixed 
From the theory, large stress concentrations 


plate. 
From the 20-percent 


appeared at the root attachment. 
chord station outboard, no large concentrations ap- 
peared. Between the root and the 20-percent chord 
station, the stress distribution was irregular so that 
coherence was not established. A smooth stress pat- 
tern was determined when these distributions were 
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Fic. 17. Stress distribution along 20-percent chord station for 


a 70-lb load applied at the apex of a partially fixed plate. Note 
Distributions outboard of the 20-percent chord station show 


similar agreement. 


plotted spanwise using the theoretical stress at the 
root, disregarding the values from the root of the 20- 
percent chord, then continuing from the 20-percent 
chord station outboard. In this manner, the stress 
pattern in Fig. 5 was produced. Fig. 15 shows the 
chordwise plot of the stresses at the 4.2-percent chord. 
Fig. 16 shows a chordwise plot of the stresses at the 
10-percent chord station compared with test results 
and the error in the theory if spanwise plotting were 
neglected. Fig. 17 shows the agreement in the correla- 
tion of theory versus test at the 20-percent chord 
station. Outboard of this station, similar conformity is 
shown. 

The spanwise plotting assumed that the theoretical 
stress distribution at the root was correct. In order to 
prove this statement, the photostress technique was 
18 displays a photograph of the 

For the particular birefringent ma- 


employed. Fig. 
method of analysis. 
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Photostress fringe pattern in the critically stressed area 
of the partially fixed plate. 
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terial used, the first fringe occurred at approximately 
6,500 psi, the second one at 17,800 psi, and the third 
one at 31,700 psi. Comparison of Fig. 18 with Fig. 5 
shows that the theoretical stress concentrations are 


reasonable. 
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Appendix A 


Rigorous Treatment of the Derivation of the 
Matrix Deflection Equation 
From the ‘Method of Analysis’ paragraph, the 
minimization of energy is given by Eq. (5): 
OU /0¢;; = —OV/0¢; (5) 


where U and V are given by Egs. (1) and (2), respec- 


vi 


7 RP tees Seer 
l = | ao(do”)*dy + f 2a;oo0" 1 ‘dy +- 5) | 2doho 2 dy + 5 } 
0 ad 0 - /J0 - /J0 


- 
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tively. In this presentation, the left-hand side of this 
equation will be presented first, and then it will be 
equated with the right-hand side. 

By taking the deflection, WW, represented by Eq. (3), 
and substituting it into the strain energy, LU, given by 
Eq. (1), the following relationship is obtained: 


I 
; l 
U : 4 } do(do")? + 2arho"di” + a2[(d1”)? + 
v4 0 


20" 2" | + 2asbi"G2” + as(h2")? + 


4ay(o2)? + 4ulacdo” + adi” + dodo” ldo + 
2(1 + u)[ao(dr’)? + 4ardr’ds’ + 4ao(po')?]{dy (15) 


7) (X) 
ar = | D(x, y)x"dx 
fo( x) 


The differentiation of the ¢’s will be done by taking 
finite differences in the following manner: 


where 


(k = 0,1, 2,3, 4) 


dysil2’ (dvi1 — On) / (16) 
on” = (dvi — 26n + bw-)/e? (17) 
Integration will be done numerically by the trapezoidal 
rule, 
N 1 VN—-1 
f f(y)dy > (f¥0 + fyy) +e pi fy; (18) 
0 2 i=1 


Applying Eqs. (16), (17), and (18) to Eq. (15) for the 
¢ terms, for example, we get 


~] 
$udydy" dod y 


= (1/2¢e%)[(1/2)ao,0(0,-1 — 260.0 + $0.1)? + do,1(0.0 — 20,1 + Go.2)? +. . -] + (1/e*)[C1/2)41,6 (G0,1 — 260,0 + 


0,1) X (di,-1 — 21,0 + O11) + A1,1(G0,0 — 260,-1 + 0,2) (G10 — 21.1 + Giz) + -. .| + (1/2¢e*) [a2,0(¢0,1 — 
2d0.0 + de) (b2,-1 — 262.0 + G21) + 2a2,1(G0.0 — 260.1 + $0,2)(G2,.0 — 2b21 + Geo) +..- | + 
(2ue*?/e*) [(1/2)ao,0(G0,-1 — 20,0 + de 1)(b2,0) + do,1(b0,0 — 260.1 + $0,2)(¢23) +...) (19) 
The boundary conditions of the problems in difference 
form are: 
OU /Odo.0 = (1/e*)[} 2a0.0(¢0,0 — Go.) + 
dy.y = O (tip singularity) } Gesldee ~ Dbes + boa)i + 
¢;0 = 0 (20) ) 2a, o(d1.0 — Ora) + A1(b1,0 — 2b1 + 1 »)f + 
ji, = YI (for j = 0,1, 2) } d2,1(b2,0 — 2¢21 + de 2) + 2a2,0(¢2,0 — $2,1)} T 


After substituting these boundary conditions into Eq. 
(19), the minimization must be performed. 


0U/d¢,,; (for j = 0, 1, 2; fori = 0,1,... 1 Vv) 


For example, 0U’/Od¢o,o is shown in Eq. (21). 


where: 


[Aool = [Do][ao][Do]” 
[Aoa] = [Do][a:][Do]” 

[Ao,2] [Do] [a2)[Do]? + 2ue?[D,] [ao] 

[Ais] = [Do]fas][Do]? + 2(1 — u)e2[D2][ao*][D2]” 


[Ai,2] 


II 


207} ao.o( —d2.0) + do1(ge 1){ |} (21) 
After minimizing with respect to 7 and 7: 
ou 1 - 1 A, 0A 9 1A, 2 \%| 
2 - [A |i); A; 0A; 1A1.2 |§¢1 (22) 
it : : AoA» 1Ao2 Is, 


[Do|[as][Do]? + 441 — u)e?[De][ai*][Do]? + 2ue?[Di\ [a] 


[Aoo)] = [Do][as][Do]? + 811 — u)e?[D2][a2*][D2]” + 4e*[@o] + 2ue*{ [D,][a2] + ({D,][a2])7} 


[A mn | [A nm - 
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where: 
p —2 1 
0 1 -—2 1 
[Do] = 1 —2 
a l 
(N) X (4) 
r—2 1 
1 -—2 1 
tale L— 4 
1 -—2 
(N) X (N) 
rl —1 
D 0 1 —1l 
[D2] = 1 1 
a | 
(N) X (N) 
do. 
0 0,2 
[do] = erage 10 
i 0 
(N) X (N) 


This completes the derivation of the left-hand side of 
Eq. (5). 

The right-hand side of Eq. (5) will now be developed 
by substituting the deflection, W, represented by Eq. 
(3), into the potential energy, V, denoted by Eq. (2). 
This yields:* 

N x 
= LD (bo + Xirbra + Xi272,)Pi2 (23) 
i=0 x=0 
Eq. (23) must be minimized for 
OV/d¢,,; (for j = 0, 1, 2) = 0,1,... ! NV) 


For example: 








rs 4 1] 0 
0 err Ff 

. | 0 

0 e| 0 

he ee, | X12 0 
0 0 X2,1 x 

[|B] = a ; | 0 

0 0; 0 

X11? V1.2 Xin" 0 

0 0 X21? Xo »2 
a 
L 0 0/0 


Eqs. (22) and (25) are substituted into Eg. (5) to obtain: 


(1/e*) [A }}os} = [B]{ Pes} (26) 
therefore, to} = e3[A]—[B]{ P} (27) 





* Subscript x refers to the number of load points along each 
station 1. 
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i 0 

0 Qk 1 

lax] = Qx,n-2 O 
g ak,N-1 

(N) x (N) 
Tax. 
. 0 ax ,2 

(a, | — Gen 10 

i 0 
(N) X (N) 

[ ax,1/2 
0 Qx,3/2 

a,*| = |- : 

la, x, nN—-1/2 





(N) X (N) 


_ fe] 


1 Py} a _ ( 
$j,Nn 
(N) X (1) 
OV/Od12 = — DL Xe2Po,2 (24) 
x=0 
Carrying out these manipulations: 
, P ; 
ol 1,z 
| |- —([B]P.2} = —[B] a > (25) 
Od; ; lp 
N,z 


(In this equation, the loading points at 7 = 0 have been 
deleted since they lie on the support and therefore have 
Where :t 


(2) 


no effect on the plate.) 


{Past = <Pis? (x) X (1) 
lp,\ 
ere... . 
1 ||. 
Ollo . 0 
0}; 1 1 0 | 
0!|0 0 
X2,7 ‘ F ¥ ‘ ‘ - | 92 NT 9 NT.. 
i...  . =! (3N) X (3Nx) 
0 | \xy.1 XN,2 ; Sew 
one... .7 
X2,2° | . . . ° . ° 
One®. :.«.8 
O | |xn,1? xy,2? et 
but {W} = [B]*{¢)} (28) 
{W} = e[B]"(A]—[B]{ P} (29) 


t Pi,o terms are omitted since these loads are on the root sup- 
port and therefore do not affect deflection or bending stress. 
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DEFLECTION AND 


Appendix B 


Rigorous Treatment of the Derivation of the 
Matrix Moment Equation 


From the ‘‘Method of Analysis’ paragraph, the de- 
flected shape of each chordwise station was given by 
Eq. (3) as: 


W = go(y) + xdily) + x%g2(y) 


The equation of the spanwise stress as given by Eq. 


(10) is: 
m Et ( 4 7) 
ee os u*) \ dy? ’ Ox? 


By taking the required derivatives of the deflection, W, 
with respect to x and y and differentiating by finite 
differences, 


Jo, a 260, i a 0, i-1 


of = 6D\(~, ¥) ) + 
2 


a 


xX 


$1,141 — 21,. + $1,i-1 


¢ 
9 Pri41 ees 2o2,; = i 2, i 1 
x" z + ude.4 (30) 
* Se f 
Inserting the boundary conditions given by Eq. (20) 
and an appropriate expression for o, (in this problem 
o, = 6M,,/t* psi/inch), Eq. (30) becomes: 


) | 
({[Ds|} dot + x?[Ds]f ait + 


[D(x, y 
{V,,:} = . 


¢ 


x?[Ds]} bof + 2ue2}¢.(MOD)}) (31) 
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where: 
2 
—2 l 
[D3] = . a 9 (N + 1) X (NY) 


\? 


}oo( MOD) g21 >(N + 1) X (1) 
feat 
Do,z 
[D(x, y)] = | 9 D,.2 


Dy,z \(N + 1) X (N + 1) 


Since D(x, y) is a function of x and y, the following in- 
formation must be given for the solution of Eq. (31). 


= D(x, y) 

0 Do,0 D;.o . Dy» 
l Do. Di, Dya 
2 Do.2 dD, 2 Dy. 


For example, D,» represents the stiffness at station 
N = lats = 3. 
inch-pounds per inch at each station 'V for every dis- 


In this manner, ./,, will be given in 
tance along N given by x. In the particular problems 
shown in this paper, the stiffness D(x, y) was assumed 
to be constant for all x at each station VV. The stress 
a, is then determined. 

A quick check on the results of the moments [Eq 
(31)] is to integrate the moment in inch-pounds per 
inch along any chord station and compare it with the 
applied moment; these should be compatible. 
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Minimum-Energy Intercepts Originating 
From a Circular Orbit 


WAYNE H. TEMPELMAN* 


Convair, a Division of General Dynamics Corporation 


Symbols 
e = eccentricity 
r = radial distance 
R = r-/ry 
v = velocity 
V = normalized velocity v/v, 


AV = normalized impulse magnitude 
normalized radial component of velocity 


a = 
8 = normalized horizontal component of velocity 
+ = angle between AV and the horizontal plane 
6 = angle between AV and the vertical plane 
= angle between circular orbit plane and transfer orbit 
plane 
0 = central range angle 
= true anomaly 
o0 = lateral range angle 
Subscripts 
¢ = circular 
f = final 
m = maximum 
o = resultant velocity vector 
p = perigee 


Abstract 


An analysis is presented of one-impulse minimum-energy in- 
tercepts which originate from a circular orbit. For the case of an 
unspecified position in the orbit, the solution involves the simul- 
taneous solution of two quartic equations with two unknowns. 
For the case of a specified position in orbit, the solution is de- 
fined by one quartic equation with one unknown. Also analyzed 
is the case of an inward intercept with a range angle in excess of 
180° under the constraint that the intercept occurs before crossing 
the desired intercept altitude. 


Introduction 


_— INTERCEPTS herein considered concern the 
transfer of a vehicle from a circular orbit to some 
point in space. The transfer is initiated by the applica- 
tion of an instantaneous rocket impulse, resulting in an 
orbital change such that the new orbit intercepts the 
desired intercept point. The minimum-energy intercept 
is that intercept which minimizes the impulse required 
to effect the transfer. 

The general problem involves intersecting a given 
point in space from an unspecified position in orbit, or, 
as the circular orbit is uniform, the general problem 
can be viewed as intersecting a specified lateral range 

Received September 26, 1960. Revised and received January 
19, 1961. 

* Flight performance and guidance analysis group, Convair As- 
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from a fixed position in orbit. The solution to the 
problem of intercepting a given point in space from a 
fixed point in orbit is obtained directly from the solution 
to the general problem. 

Assumptions are the usual two-body assumptions 
where the system consists of the attracting body consid- 
ered as a point mass and an orbiting body with negli- 
gible mass compared withthe attracting body, and where 
the only force assumed acting is gravity. 

All velocities are normalized with respect to the 
circular velocity of the circular orbit. A reciprocal 
normalized distance ratio is employed using the radial 
distance of the circular orbit. 


Statement of the Problem 


The application of an instantaneous impulse in the 
proper direction to insure intercepting the desired point 
results in an instantaneous velocity increment (AV) 
which is to be minimized. 

The variables involved in the general problem are: 
(1) a position variable (0) specifying the location on the 
circular orbit at which the impulse is applied, and (2) a 
trajectory variable (8) to represent the family of trajec- 
tories that can be passed through any two points in 
The constants involved are a distance ratio (R) 
and a lateral range (o). The variables and constants 
are discussed below and are illustrated in Fig. 1. The 
minimum AV is defined by the standard minimum- 
maximum technique of setting the partial derivatives 
of AV with respect to the variables equal to zero. 


OAV /08 = OAV/00 = 0 (1) 


space. 


The problem of intersecting a point in space from a 
fixed position in orbit with minimum impulse differs 
only in that the position variable @ is a constant, result- 
ing in a solution involving one equation with one 


unknown. 


Derivation of the Impulse Magnitude 


The application of the impulse results in an initial 
transfer orbit velocity which may be specified by the 
angle (f) between the circular orbit plane and the 
transfer orbit plane and by the normalized horizontal 


(8) and radial (a) components of velocity. The angle f 


is a function of the central angle (@) and the lateral 
range (a) (see Fig. 1). 


f = sin—! (sin o/sin 6) (2) 
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Fic. 1. Velocity and geometry nomenclature. 


The impulse magnitude is derived from the geometry 
of Fig. 1, noting that the initial normalized circular 
velocity is equal to 1. 

AV = $1+4+ a? + 8 — 28[1 — (sin? o/sin? 6) |! ot 1/2 
(3) 

The derivation is completed by deriving a general ex- 
pression for a in terms of the constants and variables 
(see Appendix). 


a = [8(cos 6 — R) + 1 — cos 6)/8 sin 8 (4) 


where KR = Vemcoras/* pinat- 


Solution to the General Problem 
The minimum impulse is obtained using Eqs. (2), 
(3), and (4), after simultaneously solving the following 
two quartic equations for 6 and @. 


0AV/0B = 0 = B4(1 + R? — 2Rcos 6) — 
63 sin? [1 — (sin? o/sin? 6)!/2] — (1 — cos @)* (5a) 





“sr 


9 R<l os cos") R 
1/2 
| R>1 @s= cos~*(#) 











180 7s 
fe) ' ov 90 


Fic. 2. Approximate locus of optimum @’s. 


OAV /06 = 0 = B (cos 8 — R)(Reos 6 — 1) — 
(8° sin @ cos @ sin? a)/(sin? @ — sin? a)'/?] — 


6?(1 + R)(1 — cos 6)? + (1 — cos 6)? (5b) 


The angles at which the impulse is applied are found 
from the geometry of Fig. 1. 

y = cos~![8? sin? f + (1 — B cos f)?]'/?/AV] 

(6) 

6 = cos~'[a? + (1 — B cos f)?]'/?/AV \ 

The nature of the equations is sufficiently compli- 
cated to warrent the use of a computer to obtain specific 
solutions. There are several restricted explicit solu- 
tions, however, which are readily verified. 

(1) Lateral range of zero degrees:* 








Fic. 3. Illustration of @ = 90° intercepts 


a= °| 
B = (2/1+R))2 AV: (4) 
\(2/(. + R)]}¥2 - || 


@ = 180 f=0 


(2) Lateral range of 90 degrees: 


f = 90° ) 


l (1+ R?2)¥2— R 


p ae - . a= : | 
(1 + R*)*/* (1+ R*)”4# > (8) 


AV = eae + R*)!/2(] 28) 2 
iad (1 + R2)"2 


6 = 90° 


(3) Central angle of 90 degrees, KR > 1: 
f = ¢ = cos 1/R'’ 6B = 1/R"? 
» (9) 
a =0 AV = [1 — (1/R)}!/*) 
(4) Central angle of 90 degrees, R < 1: 


*As A V is not a function of 8 when f = 0° and @ = 180°, 
this solution is obtained by employing the variables a and @. 
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B=1 } 
1—R AV = (3 — 2R — R%)"/\ 


f = ¢ = cos'R? 
(10) 


Qa 


For a given & ratio, the absolute minimum impulse is 
required for a lateral range of zero degrees with a down 
range of 180°. As o increases, the optimum value of @ 
rapidly approaches 90°, as illustrated in Fig. 2. The 
optimum value of @ remains slightly greater than 90° 
as o increases until the value given in Eqs. (9) and (10) 
is attained. Greater values of o result in optimum 6’s 
slightly less than 90°. When a is 90°, @ is also 90° (the 
only possible value for @). 

In general, a 0 of 90° is a good approximation to the 
optimum 6. The reason why the optimum value does 
not always occur at 90° (which is suggested by the fact 
that a 6 of 90° corresponds to a minimum plane change 
intercept) is apparent by comparing the intercept that 
occurs at 6 = 90° with one that requires aslightly greater 
plane change (see Fig. 3). The change in impulse can 
be separated into two parts, one to increase the amount 
of plane change, and the other to compensate for the 
change in range. As the slope 00/0f is infinite at 9 = 
90° (constant o), the change in impulse due to the 
change in range is predominant over the impulse change 
to rotate the transfer plane. This has the tendency to 
move the minimum energy intercept point a small 
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Fic. 5. Minimum-energy intercepts for R = 0.842 and @ = 90°. 


amount in the direction of decreasing impulse require- 
ment, determined by the sign of the slope OAJ’/08 at 
9 = 90°. Examples are shown in Figs. 4 and 5. 


Solution to the Restricted Problem 


For the case of a specified location on the circular 
orbit, the constants @ and f will be employed, with the 
variable being 8. Any point can be specified with either 
f < 90° or f 2 90°. The latter way forces the initial 
transfer trajectory velocity to have a velocity compo- 
nent along the circular orbit in the opposite direction 
to the initial circular velocity. By examining Fig. | 
and Eqs. (2) and (3), it is apparent that the same result 
could be obtained by assuming a negative value for 8 
with / restricted to £90°. The latter method will be 
employed here. 

The solution is thus obtained from Eq. (5a), after ¢ 
is eliminated using Eq. (2). This quartic equation has 
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one positive real root and one negative real root, ac- “[ | | os ee a 
cording to Descarte’s Rule of Signs. The negative real 
root corresponds to the optimum trajectory under the ad | | _ ae re 
conditions of a forced change in the orbital rotational | | | | 
direction. 
It is apparent from examining Eqs. (4) and (5) that } { } i sis | 
a 6 greater than 180° always requires the same impulse 
as for a range of 27 — 8@, the only difference being that 
a aud y change signs. T —< : ae ia, 
The restricted explicit solutions are those for f values 
of zero and 90° which are identical to the solutions given L 
in Eqs. (7) and (8) and the following solutions [exten- ° a woe 
sions of solutions given in Eqs. (9) and (10) ]. S$ 
(1) When R > 1: — 
7 =") | 
f = cos 8 = 
R — cos 6 — = 
1 — cos @ \!/? 
(i — cos ,) ew 
es a = 0 AV = | 
[1 — (1 — cos 6)/(R — cos 6) ]'/2 
"7 30 60 90 120 150 7 
—_l 9 (deg.) 
Fic. 8 Angle of impulse application for R > 1 with f/ 0° 
90°. 
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Fic. 7. Velocity required for R > 1 with f = 0°. 
(2) When Rk < 1: 
f = cos [R? — cos*@ + 
2 cos 6(1 — R)]/sin? 0 
ar (12) m 
a = (1 — R)/siné B=1 0 30 60 90 120 150 180 
4 (deg.) 
AV = {2 + [(1 — R)/sin 6]? — 2cos f} ai Fic. 9. Velocity required for R < 1 with f = 0°. 
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For R > 1 with f = cos~!{2/(1 + R)]'/*, a down- 
ward-directed application of velocity (—a@) is always 
used to obtain a minimum-energy intercept. As / in- 
creases, the range of @’s where a downward-applied 
velocity is utilized decreases (see Fig. 6) until, for f = 
90°, only upward applications of velocity are applied. 
For the latter case, a, 8, Vo, and y» are equal to those re- 
quired for a minimum-energy intercept starting from the 
initial point with zero velocity. 

For f = 0° and R > 1 (see Figs. 7 and 8) y = 0° for 
9 = 0°, increases to some maximum value less than 90°, 
and then decreases back to 0° at 6 = 180°. Increasing 
R decreases the maximum y which occurs at a greater 
range. The velocity decreases from 1 for 6 = 0° to 
1 — [2/11 + R)]!/* at 6 = 180°. 

For f = 0° and R < | (see Figs. 9 and 10) y increases 


from tan~![2(1 — R)]'/* to 180° as @ increases from 0 
to 180°. For 6 = O°, y increases with increasing R. 


As 6 increases, this situation gradually reverses, until 
for 6 > 20°, y decreases with decreasing RK. The 
velocity decreases from (3 — 2R)!/? to [2/(1 + R)]'/* 


— 1 as 9 goes from 0 to 180°. The actual transfer 
velocity (V,) exceeds the escape velocity in some cases 
(see Fig. 11). For example, V, = 1.56 for R = O and 
p= 98. 


An Altitude Constraint on the Restriction 
Problem ( R > 1) 


For 6's in excess of 180°, the previous analysis results 
in transfer trajectories which pass through the final al- 
titude before making the actual intercept. This con- 
dition is not tolerable in some cases, an example being 
descent to the surface of the earth. For the case of when 
the original orbital rotational direction is maintained, 
the transfer trajectory will remain above the intercept 
distance before making the desired intercept if 27 — @ is 
equal to or greater than the perigee angle [Eq. (A3); 
see also Fig. 12]. The constraint is 


27 — 0 


IV 


A 


p 


or cos #, 2 cos 6 (13) 


Minimizing this expression for 6 using Eqs. (A2), (A3), 
and (+) 


8 = [(1 — cos 6)/(1 — R cos 6) |!” (14) 


This value cf 6 defines the minimum energy inter- 
cept, as it results in a positive OAI’/08 slope [Eq. (5a) 
had only one real positive root]. This corresponds to 
a transfer trajectory which has its perigee located at 
the intercept point. 

An additional constraint on the transfer trajectory in 
this case is that the energy of the transfer orbit must be 
less than parabolic to avoid escape to infinity before 
making the desired intercept. Settinge = land @ = 27 
in Eq. (A2) results in 


B = (2/R)'” (15) 


Using Eq. (13) 


6< 2m — cos™'[(2/R) — 1] = 8, (16) 


The angle 6,, defines the maximum reachable @ when 
maintaining the orbital rotational direction—i.e., a 
positive component to projection of Al’ along IV, 
Whenever @ is greater than 6@,,, the orbital rotational 
direction must be reversed. As discussed previously, 
the minimum-energy solution in this case is defined by 
the negative root of Eq. (5a). 

For a given & ratio, the line separating the two re- 
gions follows the trend indicated in Fig. 13. For f = 0°, 
6, serves as the dividing line separating the reverse in- 
tercept region from the perigee intercept region. For 
6 = 6,,, the velocity requirement for the reverse inter- 
cept is greater than the perigee intercept. Increasing 


f while holding 6 = @,, results in an increasing AV for 


the perigee intercepts with decreasing Al’ for the re- 
verse intercepts. At some value of f(/,,) less than 90° 
the two AV’s are equal. Greater values of f result in 
some 0 < 6,, at which the AV’ requirements for the two 
types of intercepts are equal. 


Appendix 


Simplified Two-body Equations of Motion 


The simplified two-body system consists of an at- 
tracting body of point mass and an orbiting body with 
negligible mass compared with the attracting body. 
The only force assumed acting is gravity. An orbit is 
determined by specifying a distance (7), a velocity 
vector (V, y) and the quantity G./, which may be 
eliminated by normalizing the velocities with respect to 
the circular velocity at the initial point. The normal- 
ized velocity may be resolved into a radial (@) and a 
horizontal (8) component. 


8 = Vecosy a = Vsin y/ 


(Al) 
V =0/v, = v(r/GM)!!?_—s\ 
The general equation of the conic section is 
R = (1 + e cos ¢)/B? = Finitia/Tsinat | 
(A2) 
where e = [a’8? + (6? — 1)?]'/*\ 


The true anomaly for the initial point (referred 
to as the perigee angle) is given by (Fig. 14) 
esin 6, = a8 ecos 0, = B* — 1 (A3) 
For the purposes of intercept problems, a general ex- 
pression involving one variable is required to represent 
the family of trajectories that can be passed through 
any two points in space. The variable selected is 8, de- 
fined as positive when its component along the circular 
orbit is in the initial circular orbital rotational direction. 
A general expression for a may be derived using Eqs. 
(A2) and (A383). 


a = [82(cos @ — R) + 1 — cos 6|/8 sin 6 (A4) 








Vibrations of Honeycomb Sandwich Cylinders 


HU-NAN CHU* 
The Martin Company 


Summary 


A Timoshenko-type theory of vibrations of circular cylindrical 
sandwich shells with honeycomb cores is presented in this paper. 
The free vibrations and the propagation of free harmonic waves 
in the axial direction of an infinitely long shell are studied. 
Numerical results are given for shells with radius-to-thickness 
ratios ranging from 10 to 120 and discussed in the light of pre- 
vious works on homogeneous cylinders. 


Symbols 


x, 0, 2 = axial, circumferential, and thickness coordi- 
nates, respectively (see Fig. 1); also indicate 
quantities along these directions when used as 
subscripts 


i = indicates the layers 1, 2, 3 (Fig. 1) when used as 
a subscript 

Uri, Ugi, Uzi = axial, circumferential, and radial displacement 
components, respectively (¢ = 1, 2, 3) 

u,v, W = axial, circumferential, and radial displacements 
of shell-middle surface, respectively 

hy, he = half thickness of core and thickness of one fac- 
ing, respectively 

h =h+h 

R = mean radius of shell 

e = strain components when used with triple sub- 
scripts [first appear in Eqs. (2)] 

o = stress components when used with triple sub- 
scripts [first appear in Eqs. (6)] 

Nzi, Noi, 

Nori, Mri, } = shell-stress components defined by Eqs. (6) 

Mg, etc. | 

ri, Wr2, | = shear angles of the core and facings in axial di- 

Voi, Wor f rection and circumferential direction, respec- 
tively 

Tolacls = x, 6, and z components of .the external traction, 
respectively 

F,, Fo, g | = resultants of external tractions defined by Eqs. 

Mz, me f (7) 

FE», G2, v» = Young’s modulus, shear modulus, and Poisson’s 
ratio, respectively, of the facings 

Ki, Ke = shear coefficients introduced in equations such 
as (10) 

Grn, Gi, = transverse shear moduli of the honeycomb core 


in axial and circumferential directions, re- 
spectively 

FE, Gi, 1 = Young’s modulus, shear modulus, and Poisson’s 
ratio, respectively, of the isotropic core 





Ca, Ca = VGra/p, VGwi/pi, axial, and circumferential 
phase velocities of shear waves in honeycomb 
core, respectively 


Pi, p2 = mass density for core and facings, respectively 
Bi, Bi = w/Cu, o/Ca, respectively 

Yp» rn»\ = p2/pi, h2/h, E2/Gra, Gigi /Gra, respectively 

'r, 'G 


Received December 1, 1960. Revised and received January 
26, 1961. 
* Structures Research Staff, Denver Division. 


w, L = circular frequency and wavelength, respectively 
é = 27/L, wave number 
w/t = C, phase velocity of harmonic waves 
Q = w/t 
= th 
n = number of circumferential waves 


(1) Introduction 


F  piosegee much work has been done on the vibra- 
tions of cylinders, including composite and, in 
particular, isotropic-cored sandwich cylinders, the vi- 
brations of honeycomb-cored sandwich cylinders have 
not yet been investigated. There is a certain charac- 
teristic connected with the sandwich construction that, 
as will be seen later, has some implication for the accu- 
racy of the shell equations. By accuracy we mean that 
they are compared with the three-dimensional equations 
of elasticity; because in this respect, all shell theories 
are, of course, approximate. In order that a clearer 
account of the underlying thought can be given, how- 
ever, we shall first recapitulate briefly recent develop- 
ments of those theories of vibrations of homogeneous 
plates and cylinders that have direct bearings on what 
we shall presently discuss. 

When we say comparing the approximate equations 
with the three-dimensional exact equations, we must 
first decide on a basis of comparison. For this, Mind- 
lin' points out that in problems of small vibrations, the 
relation between frequency and wavelength is of prime 
importance. He studied this relation for the elastic 
plate in great detail by plotting the Rayleigh solution’ 
which is then used as a reference in constructing his 
approximate plate theories. In each order of approxi- 
mation, a uniqueness theorem for solutions is estab- 
lished and a self-sufficient mathematical system is pre- 
sented. His zero-order theory is equivalent to the 
classical theory of extensional vibrations of Cauchy* 
and Poisson,‘ while his first-order theory corresponds to 
the two-dimensional analog of the Timoshenko beam 
theory. Timoshenko-type plate theory was previously 
given by Uflyand.’ But the Mindlin theory is superior 
to the Uflyand theory in at least two aspects. First, 
Uflyand did not obtain the correct boundary conditions 
resulting from a uniqueness theorem, and thus found 
only an invalid solution to his problem. Second, 
Mindlin introduces a shear coefficient, x, which is de- 
termined by requiring the lowest antisymmetric thick- 
ness-shear frequency at infinite wavelength of the exact 
theory to match that of the approximate first-order 
theory. Thus, in this simple way, the approximate 
and the exact theories are beautifully connected. 
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When the Poisson’s ratio is greater than 1/3, the 
matching referred to above becomes poor in the neigh- 
borhood of the infinite wavelength. That is, the fre- 
quency curves from the approximate and exact theories 
have greatly different slopes and curvatures at the 
matching point. If we observe the plot of the exact 
Rayleigh solution in reference 1, we find that the slope 
of the lowest symmetric thickness-shear frequency at 
infinite wavelength jumps from —4/zr* to +4/7? as the 
Poisson’s ratio increases from less than 1/3 to greater 
than 1/3. In the meantime, the lowest symmetric 
thickness-shear frequency at infinite wavelength, inde- 
pendent of Poisson’s ratio, decreases from greater than 
to less than the lowest symmetric thickness-stretch 
frequency at infinite wavelength which is dependent on 
the Poisson's ratio. Thus, we might emphasize the 
fact that both these phenomena occur at the same time. 
That is, the sudden change of slope of the lowest sym- 
metric thickness-shear frequency curve and the coinci- 
dence of the frequencies of lowest symmetric thickness- 
shear and thickness-stretch modes occur at the Poisson's 
ratio, 1/3. In the second-order theory of Mindlin and 
Medick,® more «x’s were introduced to achieve a better 
matching of the slopes and curvatures. 

A similar development was carried out for the case 
of the solid circular cylinder,”* except that in this case 
the Pochhammer solution’ takes the place of the Ray- 
leigh solution. 

For the hollow circular cylinder, or the circular cylin- 
drical shell, the three-dimensional exact solution was 
obtained by Gazis" using a high-speed digital computer. 
A somewhat different version was given independently 
by Greenspon.'! Such plots, therefore, should be useful 
in constructing approximate shell theories in the same 
way as the Rayleigh solution is useful to the approxi- 
mate plate theories. Such a theory (corresponding to 
Mindlin’s first-order plate theory) had in essence been 
given by Mirsky and Herrmann,'* although phase 
velocities were plotted therein instead of the frequen- 
cies. Gazis compared his numerical results with those 
of the Mirsky-Herrmann theory and found that the 
approximate theory becomes appreciably inaccurate 
when the wavelength is of the order of the shell thick- 
ness. Gazis suggested more «’s for better matching, 
apparently in analogy with the second-order plate 
theory of Mindlin and Medick,* or the second-order rod 
theory of Mindlin and McNiven.‘ 

In order that approximate theories for the sandwich 




















zy 
a 
Padi 2 : 
Pr : ho wer 
ae 1 | 
| ih b ——“ 
i 
T 
: ho * 
Fic. 1. Cross section of shell element with sandwich-type wall 


construction. 





931 


uw 


{2 a 
g . 
oO 
iS 
» 4 
ie 
>S 4 
=| 
SI 
$|5 
S\2 , 
o|” 
F | 
Oo} 
£;°o 
a) 3 
a= 
oj-= .- 
x6 2F 
al 
> 
|@ 
3 
lo 2 
la 
re 
0.02 0.04 0.06 0.08 0.10 0.12 0.14 
5 = w Dhell Thickness 
Wave Length 
Fic. 2. Phase velocity vs. wavelength, » = 0, R/2h = 30 


cylinder can be constructed, an exact solution corre- 
sponding to that of Gazis must first be obtained. There 
would be an 1S X 18 determinant for the characteristic 
equation for the sandwich cylinder as against the 6 X 6 
determinant in the Gazis paper for the homogeneous 
shell. In addition, there are more physical variables 
for the sandwich cylinder (and more so if the core is 
orthotropic) as against the three physical variables (the 
radius-to-thickness ratio the two elastic wave 
speeds) for the homogeneous shell. Therefore, it is 
rather impractical to solve the exact three-dimensional 
equations of motion for the sandwich cylinder, even 
with the help of a high-speed digital computer. 

In view of this practical difficulty, we then try the 
next best approach—.e., the of the three- 
dimensional equations for the sandwich cylinder, for 
all wavelengths, but not for the infinite wavelength 
the corresponding frequencies of the 
the exact theories must themselves 
This was done by the author 


and 


solution 


only, since 


approximate and 
thereat. 
note. !* 

Poisson’s ratio at 
thickness-shear 


be matched 
It was also found, in reference 
coincidence of the 
thickness- 


recent 
13, that the 
lowest symmetric 
stretch frequencies was raised for sandwiches used in 
For the honeycomb-type 


in a 
and 


lightweight construction. 
core, the Poisson’s ratio associated with the thickness 
direction is zero; hence the meaningful parameter 
would be the ratio of the elastic modulus in the thick- 
ness direction to the transverse shear modulus. But it 
can easily be shown that this postponement of coinci- 
dence of frequencies also holds true for the honeycomb 
sandwich although it was not specifically mentioned 
in reference 13. In conjunction with the observations 
connected with the Rayleigh solution mentioned above, 
it is reasonable to hope that a Timoshenko-type theory 
for the sandwich cylinder would be more accurate than 
a similar type theory for the homogeneous cylinder 

pending the more exact re- 
It is further hoped that 


e.g., Mirsky-Herrmann 
sults of future investigations. 
the study of high frequency vibrations of the honey- 
comb sandwich cylinder can be facilitated by the theory 
advanced in this paper. This may have direct implica- 
tion in the design of some flight structures because one 
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possible way of improving the reliability of instruments 5 5 
lies in the understanding of high-frequency vibrations \5 ci 
of their supporting structures, which may be of the £4 e) 
honeycomb sandwich construction. 13 . 
The theory was developed by employing the varia- a > 
tional principle. The details of the derivation are =|23 
omitted since these are rather straightforward and can be Bis ; 
referred to in a recent report,'* from which this paper als 2 . : 0 
was prepared. The free vibrations and the propaga- <s ‘ 
tion of free harmonic waves in the axial direction were 2 
then studied. Numerical results are given for shells = | 
with radius-to-thickness ratios 120, 30, and 10 and dis- ) 
cussed in the light of previous works on homogeneous 0.02 004 006 0.08 0.10 0.12 014 C 
cylinders. 5 = w Shell Thickness 
Wave Length C 
Fic. 4. Phase velocity vs. wavelength, n = 2, R/2h = 30. 
(2) Shell-Stress Equations of Motion ; 
¢ 
The shell-stress equations of motion describe rela- , 
ous cylinders. 
The strain-displacement relations were also given'? as 
3 . ie Ones ait l (1. 4 “ey, 
5 ox R+2 06 ( 
; F ; ena Ou-; 
le a 
3/8 3 Cri = (OUu,;/Ox) + (O1,;/02), eri = 
3 Ou6; Lue 
siz > 3 i * Bes 06’ 
=| 8 
"E C26; = (Oue;/Oz) + 1/(R + 2) [(Ou.;,/00) — itoi| | 
3 2 
= * Hence by using assumptions (1) in Eqs. (2), we obtain: 
a l | t Crn = (Ou/Ox) + 2(OW,;/0z) 
002 004 006 008 O10 O12 014 — Cro = (Ou/Ox) — Iy(0/Ox)(a — Ye) + 
3+ 7 Nave Length 2(OW.2/Ox) (3) 
Fic. 3. Phase velocity vs. wavelength, n = 1, R/2h = 30. €rrs = OUu/OX + y(0/Ox) (Ya — Yr2) + 
2(OW»/Ox), etc. 
tions among the shell-stress components (see the list of ~ 
symbols), the inertia, and the external forces. These ry 
equations, together with the associated initial and 3 
boundary conditions, were derived!‘ by using the prin- . a 
ciple of minimum potential energy. The procedure is 3 
outlined below. alt, 
Linear distribution of displacements across the thick- s : - J Pa 
ness of each layer is assumed. Hence we have Ra = 
Un(x,0,2,t) = u(x,0,t) + sn(x,6,0) = 5 J 
Uz2(x,0,2,t) = u(x,0,t) — byWalx,0,t) + 5 06 - 
( + hi) P22(x,8,0) é 2 
U73(x,0,2,t) = u(x,0,t) + Iivalx,0,t) + « 
(f — hy) W22(x,6,t) “a o* 
e(x,0,2,t) = w(x,0,t) (1) = 
Uo (x,0,2,t) = v(x,0,t1) + sho(x,6,0) % 3 
Ue2(x,0,2,t) = v(x,0,t) — hy er(x,0,t) + we 02 i 
(s + hy) Woo(x,0,t) = . 4 
U63(X,0,2,t) = v(x,0,t) + hyer(x,0,t) + a | 


(s — hy) Woo(x,0,t) 01 02 03 04 05 
5 = Shell Thickness 


which may be regarded as a natural extension of those Wave Length 
2, R/2h = 120. 


assumed by Mirsky and Herrmann’? for the homogene- Fic. 5. Frequency vs. wavelength, n = 0, 1, 2, 
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With the displacements (1) and the strains (3), we 
can calculate the strain energies, kinetic energies, and 
external work, and can perform the variation similar 


ON, 1 ON¢- 
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to that of reference 12. 
The shell-stress equations of motion are then ob- 


tained as 


: h;? . hh» ; 2 
+ F, = 2phi\ ud + vn) + po} 2h + Rp et In) a — Yn) + oF (h? — In*)bz2 
\ Ol 


») 


(h® — hy ve. 


Oo’, 1 OMe, 2 2 it 2... ul 
+- — Or + mz = = phi®| dn + R + po} Asholh + n)(¥n — vez) + (Ch? —In®)\ Ye + R 
o \ o \ 


Ox R 06 


OM, fe) : ‘ 1 | OMe, 
+ hy ( N i N 7) + | + hy ( N 623 ~~ 
R O86 


Ox Ov 06 


2 j (v 4 ut 
Pilly” | Yr : 
3 ei 


OM, 1 Oo. 2 d 
6 6 _ Oo + me = os (Uo + ) + pa| Ia + hy) (va 


Ox R o6 ; R 


06 O00 Ox 


R 


where V, = > N,i, ete. (5) 
i 1 


and 

h “ 

Na = a ie —) dz 

; i h __ ( + *) ‘ 

- hy oe 

Nw = | zr2\ 1 = Iz j 
= : ( i) (6) 
h ™ 

Nz = [ Orr3 (: + ae etc. 


are the shell-stress components. The o’s in the above 
integrands are the stress components one-to-one corre- 
sponding to the strain components of Eqs. (3) (identi- 
fied by their subscripts). 

In Eqs. (4), F;, Fe, mz, me and qg are the external 


forces, and are defined by 


iS z=h 
F,, Fg, , = fe fo, aC + all = 
eee (7) 
m:,m = E fez (1 ar “) | | 


The initial conditions are the specifications of u, v, w, 
Yn, V2, Wer, Wo2, and their first time derivatives through- 
out the shell. 

For boundary conditions, we have: 

(1) Along a longitudinal boundary @ = constant, one 
member of each of the seven products Neu, Nev, Qew, 


Nor2) |- On ~ 


) + pf Bahn one 


1] OM re) OM, re) . ae 
| _ + hy (Nes — v0) | + > + hs 5 (N63 — N02) — Va + 
: x 


2 3 
e piti*( Un + ) + po Bex 
2) R 


(4) 


hy 
Mm, = 
h 


V2) + hyheo(h + hy(v, +") 
R 


2 v 
— We) abe (h? — in)( Ue T )| 
3 R 


hy ( QO O ) m hy 
oo 2 m = 
oe = cs 


— V2) + Iyho(h oa hy) (ve + *) 
\ 


Mowe, (Mon + Ii(Now — Nor) (Wn — vr), Mever, 
[Mo + hi(Nes — Nez) \(ve — We) must be specified. 
(2) Along a circumferential boundary x = constant, 
one member of each of the seven products V,u, N,v, 
O,w, Mx, [Mi + hi(Nz3 — Nw») \(vn — Wr), M ,oe2, 
[Mero + Ai( Nis — Naor) |\(wor — Ye2) must be specified. 
These initial and boundary conditions may be shown 
to be sufficient for the uniqueness of the solution follow- 
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ing reasoning similar to that of references 8 and 12. 
When h, — 0—i.e., for vanishing facings—the fourth 
and fifth equations of set (4) become identical to each 
other, the sixth and seventh equations of set (4) be- 
come identical to each other, and set (4) reduces to that 
derived by Mirsky and Herrmann!’ for homogeneous 


cylinders. 


(3) Stress-Strain-Displacement Relations 


We consider only the orthotropy in which the axes of 
elastic symmetry coincide with the coordinate axes 
used in this investigation, namely, the x, 0, and z axes. 

The stress-strain relations for the orthotropic core 


may be written as 


Fa : 
Orn = (€rr1 + Ver 061) + 
1 — verre 
Ex Ver — Vor 261 
. zz1 
Ea I VoriVr61 
Ea \ (8) 
ay = = (Coo + Vrmerr) + 
1 — vornvrer 
Fea vie + veri¥r61 
. z2l 
Ey a VoriV rl 
orn = Groiera, zx = Graera, O21 = Gzme261 | 


in which e.., has been eliminated algebraically among the 
six stress-strain relations required for a three-dimen- 
sional solid. 
For the two isotropic, identical facings, the stress- 
strain relations are simply 
Fo Vo 
~ (Cxzi = V2 901) + | a 


= ’ Cust; CL. 
1— V2~ 


Crrt = 
V2 


(¢ = 2,3) (9) 


Next, the strain components of Eqs. (3) are used in 


SCIENCES—DECEMBER 1961 
qs. (8) and (9). Then the resulting stress-displace- 
ment relations are integrated through the thickness ac- 
cording to Eqs. (6) to obtain the shell-stress displace- 
ment relations. In the process, the integrals contain- 
ing o,,; are neglected and four x«’s are used in the trans- 
verse shear terms for exactly the same reasons as in 
references | and 12. 

There are 30 shell-stress displacement relations, as 
there are 30 shell-stress components of Eqs. (6). A 
typical shell-stress displacement relation is 


ne Ow — Ow ; 
Q, = KVGra2hy + Vu + K3"Ge2he + Wr (10) 
Ox Ox 


The complete set of shell-stress displacement relations 
may be found in reference 14. In this set, the approxi- 


mations 


log [((R + hy)/(R — ky)) = 2h4,/R[1 + (A?/3R?) | 


R—-—h R+h 2h» 2(h® — h,*) 
log ] ~ - 
" (j =f ) + log € + ~) R + 3R3 


fe a 2 fs + 4 h? — hy? 
og — = 
"la 5 ENR +h R2 


have been used. 


(4) Shell-Displacement Equations of Motion for 
Sandwiches Used in Lightweight Construction 


In sandwiches used in lightweight construction, the 
facings are often thin enough in comparison with the 
core to permit neglect of their flexural rigidity, trans- 
verse shear, and rotatory inertia.')'® Then Eqs. (4) 
and the shell-stress displacement components such as 
Eq. (10) can be simplified. If we substitute these 
simplified shell-stress displacement components into 
the simplified version of Eqs. (4), then we get the fol- 
lowing set of shell-displacement equations of motion: 
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describe an orthotropic cylinder. If further the 
orthotropy tends to the isotropy, then Eqs. (11) be- 
come exactly those derived by Mirsky and Herrmann'* 
for the homogeneous isotropic cylinders. 

The Mirsky-Herrmann theory has been shown by 
Gazis'® to be valid for shells as thick as one-quarter of 
the mean shell radius if the wavelengths are long in 
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comparison with the shell thickness. Gazis also showed 
that for shorter wavelengths (down to about the shell 
thickness), the fourth and fifth modes (for » < 0.3, 
identified as longitudinal and circumferential shear, 
respectively) deviate appreciably from the exact theory. 
He attributes this deviation to the coupling effect from 
the higher thickness-modes. In a second-order plate 
theory, Mindlin and Medick® had pointed out that 
coupling from the higher thickness-modes becomes ex- 
tensive when the Poisson's ratio increases to a value at 
which the first symmetric simple thickness-shear and 
stretch frequencies become coincident. 

In sandwiches used in lightweight construction, this 
Poisson’s ratio at coincidence is raised.'* Thus, if the 
frequency curves from the three-dimensional exact 
theory for sandwich cylinders were plotted, we could 


expect that the changes of slopes and curvatures at 
infinite wavelength would not be as rapid as that ob- 
served in the plot of the Rayleigh solution.! Hence, 
it is reasonable to hope that the theory developed in 
this paper would be more accurate than the correspond- 
ing theory for the homogeneous cylinder,'? pending the 
results of further investigation. It is believed, how- 
ever, that frequency curves corresponding to that for 
the Rayleigh solution (reference 1, Fig. 2.114) should 
first be obtained for the sandwich plates. This is now 
being contemplated for a forthcoming report. 

Gazis’ results were chosen for comparison rather than 
the other three-dimensional solutions'' because the ap- 
proach used in this paper makes a comparison with 
Gazis more convenient. A comparative study of vari- 
ous approximate and exact theories for the homogene- 
ous cylinder has been made by Greenspon." 

When (1/R)? < 1, Eqs. (11) can be further simplified 
for both isotropic-cored and honeycomb-cored sand- 
wiches. In honeycomb sandwiches used in lightweight 
construction, E,,//2 and E/E: are usually of the order 
of 10~‘ to 10~, while /,/h2 is usually of the order of 10. 
Therefore, terms containing E,, Ke, and G,» can be 


dropped. Hence, we have 
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Eqs. (12) break up into two uncoupled sets when 
O( )/06 = 0. One set of three equations describes 
the axisymmetric motions while the other set of two 
equations describes the purely torsional motions. 
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(5) Solution for Free Vibrations of the Infinite 
Shell for Sandwiches with Honeycomb Cores 


Solution examples will now be given for the honey- 
comb sandwich cylinder. 

For free vibrations of the infinite shell, we set F,, Fs, 
g, m,, and mg of Eqs. (12) to zero and assume solutions 


of Eqs. (12) in the form of 


u(x, 0, t) Ue —™ cos né 
Yul(x, 6, t) = We" © cos nb 
v(x, 6, t) = Ve sin no (13) 


Yal(x, 0, tf) = Wee ©) sin 0 
w(x, 0,t) = Weis! — 


where 7 = 0, 1, 2, ... is an integer indicating the num- 
The following charac- 


cos n@ 


ber of circumferential waves. 
teristic equation is obtained: 


ay Q12 13 a4 15 | 

ay 22 23 24 25 

13 23 33 34 a3, = O (14) 
ay4 4 34 44 45 

| dis 5 35 45 155 | 


where 





ay) 


W 





ite 
res 
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2E she 2Goh: 
g 


a = 2(pily oa pohts)w = ri 2 = n*, etc. 
1 — pw°* R? 
(15) 

When n = 0, Eq. (14) reduces to 
ay” ay” a5" 

433° Q34°° 
ay,” Qo" 5° = 0 (16) 

34° 144 0 
a5" 2% - 155 . 
where ay” = lay = etc. 


The 3 X 3 determinant of Eq. (16) describes the 
axisymmetric motions. The 2 X 2 determinant of 
Eq. (16) describes purely torsional motions that are 
now uncoupled from the axial and radial motions. 
Further, when & = O0—i.e., for infinite wavelength 
the 2 X 2 determinant of Eq. (16) furnishes the two 


cutoff frequencies of the torsional modes. They are 


qg= © (circumferential) 
kK. | i. (circumferential (17a) 
20 nil thickness-shear) 


MN (1 + 37,rp) 


where the /,.? term has been neglected because of the 
thin facings assumption. 
On the other hand, the 3 X 3 determinant of Eq. 


(16) for — = O gives the following cutoff frequencies: 


= @ (longitudinal) 
1] Eoho (ring extensional, 

“- RN i-tiethas @ ““breathing’’) 

e.. _ (longitudinal 


KI 
thickness-shear) 


I= | 
‘“ PA + 3 ,0n) 
(17b) 


where the /,” term has been similarly neglected 
If instead of letting m = 0, we put — = Oin Eq. (14) 
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8, R/2h = 120 


i.e., the frequency equation is solved at infinite wave- 
length—then Eq. (14) yields the five cutoff frequencies 


for the five modes of the shell theory. They are given 
by 
33.0 Asai 35 
ai Ai20 
13410 A440 450 =x @ (1S) 
A120 A220 
1350) 45.0 15510 


where di,.0) = [ay le _o,ete. The2 X 2 determinant in 
Eq. (18) describes motions corresponding to the longi- 
tudinal shear vibrations of the exact theory of Gazis, 
while the 3 X 3 determinant in Eq. (18) describes that 
corresponding to the plane-strain vibrations. The two 
kinds of motions are, therefore, uncoupled at — = 0, as 


was already noted by Gazis. 


(6) Determination of Shear Coefficients «, 
and K3 


Strictly speaking, the correct way to determine «x; and 
k, is to match the cutoff frequencies of the thickness- 
shear vibrations from the shell theory to the corre- 
sponding frequencies from the exact theory. This will 
necessarily involve the number of circumferential 
waves, m, as can be seen from Eq. (18). That is, « and 
ko must be determined for each n. Also, the exact 
theory must first be solved at infinite wavelength, in- 
cluding the effect of the number of circumferential 
waves, 1. 

This method of determining x's may prove quite 
tedious if carried out, thus destroying the original pur- 
pose of introducing these coefficients, namely, to pro- 
vide a link between the approximate theory and the 
exact theory in a simple way. Fortunately, we see 
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from reference 12 that the way n is introduced in de- 
termining the shear coefficients is through the term 
proportional to (h/R)*?. In our theory—Egs. (12) 
(h/R)* is assumed negligible compared to unity. 
Hence, consistent with our assumption, we may deter- 
mine the shear coefficients at m = 0. This is a con- 
siderable simplification. 


A. The Longitudinal Shear Coefficient *, 

We shall determine «x, by matching the cutoff fre- 
quency of the longitudinal thickness-shear vibrations 
from the shell theory given by the third equation of 
(17b) to the frequency of the first longitudinal anti- 
symmetric thickness-shear mode of the exact theory. 
From reference 13, the exact theory gives the following 
equation for the longitudinal antisymmetric thickness- 


shear nodes: 
Bihy tan Bihy = L/rpr, (19) 


Then x, can be found by equating the lowest root of 
Eq. (19) to the third equation of (17b). 


B. The Circumferential Shear Coefficient *, 

k. is determined by matching the cutoff frequency of 
the circumferential thickness-shear vibrations from the 
shell theory given by the second equation of (17a) to 
the frequency of the first circumferential antisymmet- 
ric thickness-shear mode of the exact theory. From 
reference 13, the exact theory gives the following equa- 
tion for the circumferential antisymmetric thickness- 


shear modes: 
Bi, tan By = 1/ryr, (20) 


Then «2 can be found by equating the lowest root of 
Eq. (20) to the second equation of (17a). 
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An examination of Eqs. (17a), (17b), (19), and (20) 
reveals that «; and xy will always be equal to each other. 

A remark must be added here, however. In deriving 
Eqs. (19) and (20), the thin shell assumption (/R) < 
1 had been used in reference 13. However, it can be 
shown that the equations for determining « and x2 are 
unchanged if the assumption (/R) < 1 is relaxed to 
read (h/R)? < 1. 


(7) Numerical Examples and Discussions 
These numerical data are used for the calculations: 


R/2h = 120, 30 and 10; 1, = 2/23; r, = 28.8; re = 

E>/Gin = 166.7; re = G.n/Gra = 0.6; v2 = 0.32. 
The shear coefficients were found to be 

Ky? = Ke? = 0.9975. 

The values R/2h = 30 and 10 were used in references 
10 and 12, while the ratio of the core transverse shear 
moduli, 0.6, is representative of the commonly used 
metal honeycombs. 

The calculated results are plotted in Figs. 2 through 
16. The first three figures are phase velocities vs. 
wavelengths, while the last twelve are frequencies vs. 


wavelengths. Each figure contains five curves 
(branches) since this is a five-mode theory. 
For n = 0, the torsional motions which are depicted 


as the second and fourth branches, are of course un- 
coupled from the axisymmetric motions. The plots 
of phase velocities for m = 1 and m = 2 are similar to 
that of Fig. 11 in Gazis’ paper. If we compare Figs. 3 
and 4 of this paper with Fig. 11 of Gazis’ paper (since 
R/2h = 30 in all three), we see that the maximum and 
minimum of the first branch occur at about the same 
wavelengths in both papers. Only the magnitudes are 
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different. Thus, the effect of the addition of two thin 
and stiff facings in the shell is to slightly increase the 
axial wave velocities, if the core material is used as a 
reference. This increase can be shown to depend on 
the ratio of the elastic moduli and r,r,._ Since these 
were kept fixed for all the calculations of this paper, the 
maximum phase velocities of the first branch are about 
the same for all values of R/2h. These maximums can 
be shown to vary as ¥/7,, other things being equal. 
For n > 2, all branches tend to infinity as 6 tends to 
zero, just as in the case of homogeneous cylinders 
(refs. 10 and 12). 

Figs. 5 through 16 show that for shorter wavelengths 
and higher frequencies, all corresponding branches tend 
to be the same for all » and all R/2h. The difference 
is, therefore, confined to a region near the origin. For 
n > 2, every branch has a finite frequency at 6 = 0. 
Also, as 6 increases—i.e., as the wavelength shortens 
the second and fourth branches become parallel to 
each other and the third and fifth branches become 
parallel to each other. This must be the case, because 
for shorter wavelengths, the cylinder behaves like a 
plate and all branches should tend to that describing 
the vibrations of a sandwich plate. In the plate, the 
extensional motions (corresponding to the second and 
fourth branches) are uncoupled from the flexural mo- 
tions. 

At 6 = O, the fourth and fifth branches, respectively, 
give the simple circumferential and simple longitudinal 
thickness-shear frequencies. The circumferential 
thickness-shear frequency is the lower of the two be- 
cause G.9, < G; is used in the calculations. If trans- 
verse isotropy (in which G. = G,.) were assumed, 
the fourth and fifth branches would coincide at 6 = 0. 

As to the nomenclature of the various branches of 
the frequency curves, we prefer the Gazis idea (ref. 
10)—i.e., they should be identified in the neighborhood 


of the infinite wavelength, 6 0. Thus, the first 
branches of nm = 0, 1 and m > 2 are, respectively, the 
axisymmetric shell-flexural, the beam-flexural, and the 
ring-flexural modes. The second and third branches 
are, respectively the circumferential-shear and ring- 
extensional modes. The fourth and fifth branches 
have already been discussed. 

The crossing of the third and the fourth branches for 
nonzero n at 6 = 0.296 gives an illusion that the tor- 
sional modes may be always uncoupled from the rest. 
Actually, Eq. (14) cannot be factored for nonzero n. 
The appearance of the crossing of the branches is the 
result of the choice of the numerical data. For ex- 
ample, at 6 = 0.296, Eq. (14) has two roots so close to 
each other that they cannot be distinguished on the 
scale used in the drawings. However, the region near 
the origin is clearly indicative of the coupledness or 
uncoupledness of the branches. 

Generally speaking, the overall behavior of the vari- 
ous branches in the longer wavelength region is similar 
to that of the homogeneous cylinder (refs. 10 and 12) 
This in turn gives a justification of the numerical data 
adopted--1.e., the emphasis of the variation of the 
R/2h ratio—since we are attempting to establish a 
shell theory which is expected to be valid for moderately 
thick shells. 

The curves for (2/2h) 
son) indicate the necessity of plotting the region very 
near the origin for thin shells. The case of nm = S may 
be regarded as representative of those of large n, and 
also perhaps of interest in flutter analysis when the shell 


120 (included for compari- 


is thin. 

For thin shells in which 2h/k < 1, Eqs. (12) can be 
further simplified. In thin shells, it is customary to 
measure wavelengths by the radius, 2, of the cylinder. 
In other words, only a small region near the origin is of 
the most interest. Then the transverse shear deforma- 
tion becomes insignificant under ordinary circum- 
stances and we have the commonly known three-mode 
theory for the thin shells. 
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Pressure-Gradient Effects on the Preston 
Tube in Supersonic Flow' 


J. F. NALEID* ano M. J. THOMPSON** 
United States Air Force Academy and The University of Texas, Respectively 


Summary 


This paper is concerned with an experimental investigation of 
the effects of a longitudinal pressure gradient in a supersonic 
stream of air over a bounding surface on the performance of a 
Preston or impact-pressure tube at the surface. Evidence is pre- 
sented which indicates that for the Mach number considered and 
for the range cf pressure gradients covered, the Preston tube func- 
tions in a completely satisfactory manner for the determination 


of local shear stress 


Symbols 

d = outside diameter of impact probe 

fi, fg = functions of Mach number in shear stress-probe ve- 
locity relation [Eq. (3)] 

my = Mach number parameter [see Eq. (2)] 

Mf, = Mach number outside of boundary layer 

p = static pressure in boundary layer 

/ = temperature recovery factor 

Ra = Reynolds number based on probe diameter 

Re = Reynolds number based on boundary layer momentum 
thickness at observation station 

t = ratio of inside and outside diameters of probe 

u = local velocity in boundary layer at distance y from 
boundary 

uy = free-stream velocity at outer edge of boundary layer 

Us = velocity at center of impact probe, y = d/2 

u* = friction velocity in boundary layer 

\ = distance along boundary surface 

y = distance normal to boundary surface 

+ = ratio of specific heats for gas (y = 1.400 for air) 

n = friction-distance parameter 

y = kinematic viscosity of gas 

p = mass density of gas 


a = Mach number parameter [see Eq. (2)] 
shear stress at boundary surface 


r = 
? = ratio of flow velocity to friction velocity 

?, = value of ® at outer edge of boundary layer 

w = exponent in viscosity-temperature relation (w = 0.768 


for air) 


(1) Introduction 


Des RECENT YEARS there has been considerable 
interest in the direct determination of skin-fric- 
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tion forces in turbulent boundary-layer flow over a wide 
range of free-stream speeds. This interest arises in 
connection with the basic relation for velocity within 
the boundary layer which is generally assumed to con- 
form to ‘‘the law of the wall,” originally proposed by 


Prandtl—that is, 


u/u* = f[(yu*)/v| (1) 


Here uw is the local velocity at a distance y from the 
boundary surface, «* is the friction velocity V t0/p, 
v is the kinematic viscosity of the fluid, 7) is the shear 
stress at the boundary surface, and p is the mass density 
of the fluid. In pipe flow the shear stress may be de- 
duced from measurements of pressure drop, but in the 
case of a body immersed in a fluid, this procedure is not 
possible. 

In experimental studies of this problem, at least three 
methods have been employed. The first of these is the 
“floating element” balance which employs an isolated 
portion of the area of the boundary surface, supported 
in such a way that the shear force acting on it may be 
measured directly. This method was first developed 
by Dhawan! and modified and extended by Weiler and 
Hartwig.” The latter developments have been used 
extensively at the Defense Research Laboratory, The 
University of Texas, as a basis for studies of turbulent 
boundary layers at supersonic speeds,* both on wind- 
tunnel models and on flight vehicles.‘ The flight ex- 
periments involved the use of an inertia-compensated 
balance design which cancels out the effects of longi- 
tudinal and lateral accelerations of the test vehicle. 

The second method may be described as a “‘heat trans- 
fer’’ procedure and was developed initially by Ludwieg® 
and used by Ludwieg and Tillman® in an investigation 
of turbulent shear stress at low speeds. The procedure 
involves the local measurement of heat transfer through 
an element of the boundary surface and the correlation 
of this result with shear stress by means of the Reynolds 
Analogy. 

The third method, considered here, is the use of a 
total pressure tube resting on the boundary surface, as 
originally proposed by Preston.’ The validity of the 
method for pipe flow was demonstrated by Preston and 
extended to external boundary layers by Hsu,’ both of 
these investigations being restricted to low speeds. 
An extension to supersonic flow was made by Fenter 
and Stalmach® and indicated satisfactory performance 
under these conditions when no pressure gradient existed 
in the main flow direction. 
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Schematic arrangement of DRL variable Mach number 


wind tunnel 


The limitations in the case of the shear force balance 
are mainly of an instrumentation type, although when 
used in a flow with an external pressure gradient, some 
difficulty may arise because of internal flow through the 
balance generated by pressure differences over the gap 
between the shear element and the surrounding bound- 
ary. When applied to bodies with high surface tem- 
peratures, there may be complications resulting from 
large expansion of the shear disk and binding against 
the walls of the opening in the surrounding surface. 
The mechanical elements in the skin friction balance 
also make it more susceptible to erratic behavior caused 
by vibrations and shock loadings. 

The heat-transfer method does not lend itself too 
readily to high-speed flows where heat transfer itself 
may be a part of the investigation, since the relation- 
ship between shear stress and heat-transfer rate must 
have already been established. 

The Preston tube is based on the assumption of ef- 
fectively complete immersion in the region where the 
“law of the wall’’ is applicable. It also has the ad- 
vantage of markedly simpler construction and opera- 
tion than is the case with the other instruments men- 


tioned. 


(2) Experimental Facilities 


The present investigation is essentially an extension 





DRL variable Mach number wind tunnel. 


Fic. 2. 
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of the work reported in reference 9 to include a limited 
study of the effect of pressure gradient. The DRL 
variable Mach number supersonic wind tunnel was 
modified so as to permit the development of a controlled 
pressure gradient on the floor of the test section. This 
facility was originally developed as the basis for a pilot 
study of the flow characteristics of a variable-geometry 
supersonic nozzle, involving a mechanically simple con- 
tour-changing mechanism. The basic principle in- 
volved, first proposed by Orlin,'® was to connect the 
blocks forming the throat portion of the nozzle with the 
test section by means of flexible plates. The original 
configuration is shown schematically in Fig. 1. A 
photograph of the original test section with one side wall 
removed is shown in Fig. 2. Here the downstream 
portion of the nozzle blocks and the pivot arms are 
The black circle at the left 
test section for 


clearly visible at the right. 
rear wall of the 


is an insert in the 





Fic. 3. Schematic diagram of pressure gradient test setup 


mounting pressure probes. 

A change in nozzle area ratio and Mach number is 
accomplished by rotating the nozzle blocks so as to alter 
the duct height at the throat. The flexible plates are 
unsupported between the throat and the test section 
and are permitted to bend freely as determined by 
their structural characteristics and end conditions 
The tunnel is two-dimensional, the duct width being 
The dimensions at the test section are 
The tunnel is capable of continuous 


held constant. 
2.00 in. by 2.00 in. 
operation with a battery of reciprocating-type compres- 
sors having a total power requirement of 440 hp. 

An extensive flow survey'! was conducted in the test 
section of the wind tunnel to determine the flow charac- 
teristics, both at the test section and through the 
nozzle. Reasonably uniform flow was developed for 
test section Mach numbers ranging from 1.80 to 3.80, 
the results indicating a maximum deviation in Mach 
number of 0.05. 


(3) Modification of Wind Tunnel 


The earlier tests of Stalmach, reported in reference 9, 
had utilized the wind tunnel just described without 
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Modified 2 by 2 in. variable Mach number wind tunnel 
with skin friction balance and instrumentation installed 


Fic. 4. 


major modification except for the insertion of a survey 
probe and skin friction balance in the tunnel test section 
floor. The investigations were primarily concerned with 
the effects of Mach number on the performance of a 
Preston tube mounted on the floor of the wind tunnel. 
The present study was concerned with a determination 
of the effect of a longitudinal pressure gradient on the 
tube behavior. 

The upper wall of the wind tunnel in the test section 
region was modified as shown in Fig. 3. An additional 
flexible plate was placed between the original nozzle 
wall and a drive mechanism previously used for an 
adjustable diffuser. Up-and-down motion of this drive 
produced different 
auxiliary plate and corresponding variations in pressure 
along the bottom wall or floor of the wind tunnel. 
Other details shown in Fig. 3 include the impact pres- 
sure survey system and the skin friction balance with 
its calibration device. It should be noted that both 
of these units were installed so as to provide for pressure 
probe readings at the center of the skin friction balance 
disk in the tunnel floor, although shown separately in 
Fig. 4 is a photo- 


amounts of curvature in this 


Fig. 3 for purposes of clarification. 
graph of the modified tunnel test section with the skin 
friction balance installed. 

Measurements of the static pressure along the wind 
tunnel floor under the auxiliary plate demonstrated 
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Fic. 5. Schematic arrangement of a Preston tube. 
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that constant values of the pressure gradient dp/dx 
could be developed over an appreciable length. At the 
same time a uniform distribution of Mach number across 
this region was maintained. At a test section Mach 
number of 2.00, pressure gradients ranging from dp/dx 
= —0.150 psi/in. to +0.700 psi/in. were established. 
A brief study was made as to the effect of these 
pressure gradients on the performance of the skin fric- 
tion balance, with the conclusion here that this effect 
was negligible. 


(4) Plan for Experimental Program 


After having completed the necessary tunnel modifi- 
cations to produce the desired pressure gradients, the 
actual investigation of the surface probe or Preston 
tube was initiated. As described in references 7-9, the 
Preston tube is a simple Pitot or impact tube set against 
the boundary surface and facing into the flow, as shown 





” 


Fic. 6a. Compressible flow law of the wall velocity profiles 


in Fig. 5. With the assumption that the tube is within 
the region in which “‘the law of the wall’’ applies, a read- 
ing of the total pressure at the tube opening can be 
used to determine the local shear stress at the boundary 
surface adjacent to the tube nose. In order to verify 
this hypothesis for the case of adverse pressure gra- 
dients, the following test procedure was proposed. 
After establishment of the desired values of dp/dx and 
M,, readings were to be taken of the pressure acting at 
the nose of the Preston tube. In the present installation 
the tube was mounted on a supporting strut so that it 
could be made to traverse the entire boundary layer. 
When in contact with the tunnel floor, it functioned as 
a Preston tube. 

In addition to the impact pressure measurements, a 
skin friction balance having a disk diameter of 0.25 in. 
was installed at the same location, the pressure survey 
mechanism having been removed. Independent direct 
measurements of shear force were then made under the 


same flow conditions. As in the uniform pressure tests 
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Compressible flow law of the wall velocity profiles. 


Fic. 6b 


of reference 9, the nose of the pressure probe and the 
center of the skin friction balance were located in 
identical positions. The test conditions were always 
controlled so as to hold both Reynolds number and 
Mach number at constant values during the entire test 
program. The nominal Mach number was M/,; = 2.00 
and the Reynolds number, based on momentum thick- 
ness at the observation station, was R, = 1.05 X 10# 


(5) Analysis of Results 


Following Fenter and Stalmach, the total pressure 
probe data were used to determine the local velocity 1 
within the boundary layer. These results were com- 
pared with ‘‘the law of the wall,’ writing the latter in 


the form 
(, Vo) sin" |(@/%)Yo] = f(n) (2) 


where ® = u/u*, ®; = u;/u* and wu, is the free-stream 
velocity. The parameter ¢ is a function of Mach num- 
ber of the form o = m,/(1 + m), where m, = (y — 
1)M,?/2. The quantity 7 is the friction-distance 


parameter, defined as » = (yu*)/v. The function f(y) 


«ere 








Fic. 7. Correlation of surface impact probe data with direct skin 
friction measurements for a probe diameter of 0.020 in. 


in Eq. (2) is represented by the numerical data tabu- 
lated by Coles.'* 

A comparison of experimental and theoretical veloc- 
ity distributions as shown in Figs. 6a and 6b for the 
eight different values of dp/dx. It appears that the 
agreement over a limited portion of the boundary layer 
is quite good. As in the case of the tests without 
pressure gradients, the ‘‘velocity defect law’ takes over 
when y or 7 is large. Comparisons with the latter are 
not presented here but are given in reference 13. A 
sliding vertical scale is used in Figs. 6a and 6b in order 
to facilitate the comparison for each dp/dx value con- 
sidered 

The skin friction balance measurements made possible 
the determination of the wall shear stress 7) and the 
local skin friction coefficient C, These results are 
compared with the relation given by Fenter and Stal- 


mach, which is of the form 


[fi(Mh) Ra/ Vo] sin (uo o/i1) 
V fo( My) RoC 7/2 f(V fol Mi) RC 1/8) (3) 





Comparison of test data with data from reference 9 


Fic. 8 
Here 

fi, (AN) (1 + myr)—“ te 
and 

fo( My) = (1 + myr)—“ + 
and Ra = pli; M1 


is the Reynolds number based on the outside diameter 


of the probe. Also, 
u; = velocity determined by probe reading, for 
y = 6/2. 
ry = temperature recovery factor, taken as 0.90. 


w = exponent in viscosity-temperature relation, 
taken as 0.768. 


This relation is based on the idealized limiting case 
where the inside diameter of the Preston tube ap- 
proaches zero, the wall thickness is equal to the outside 
radius, and the tube then reads the impact pressure at 


y = d/2. The effect of variations of the ratio of inside 
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to outside diameter, ¢ = d;/d,, on the probe reading 
was discussed in reference (9). Consideration of the 


extreme case of zero wall thickness, that is, ¢ = 0, 
showed no significant difference in the pressure-shear 
stress relation. 

The results of the computations based on Eq. (3) are 
shown in Fig. 7, with each point representing the re- 
duced test data for a particular value of the pressure 
gradient. The points naturally fall in a restricted range 
since only one value of 1/; = 2.00 was included in the 


investigation. 


(6) Conclusions 


A comparison with Fenter’s and Stalmach’s results® 
for dp/dx = 0 and for an extended range of Mach num- 
bers is shown in Fig. 8. It is concluded from these 
results that the probe measurements are of acceptable 
accuracy when made in the presence of moderate pres- 
sure gradients, although it is realized that the present 
range of pressure gradients is somewhat limited, par- 
ticularly in the direction of negative values of dp/dx. 
Further extensions of this study to include heat transfer 
and surface roughness effects are now in progress at 
DRL. An expansion of the program to include other 
Mach numbers and a wider range of pressure gradients 
is also contemplated. 
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The Effect of an External Supersonic Flow 
on the Vibration Characteristics of Thin 


Cylindrical Shells 


H. M. VOSS* 


Boeing Airplane Company 


Summary 


The problem is first set up, based on Reissner’s shallow- 
eylindrical-shell vibration theory and piston-theory air forces. 
rhe nature of the solution and its critical dependence on fre- 
quency are illustrated by a simple two-mode solution. The com- 
plexity of the frequency spectrum previously noted by several 
iuthors is emphasized in the inability to choose modes of maxi- 
mum coupling as is commonly done in other aeroelastic problems 
Based on these results, the problem is recast in terms of the zero- 
iir-speed frequencies. The solution is then outlined in terms of 
Gol’denveizer’s equations for thin cylindrical shells, including 
mid-plane inertia, and determination of frequencies and damping 
f the modes as a function of the aerodynamic parameters. An 
example is included which illustrates the above difficulties and 


shows the effect of structural and aerodynamic damping 


Symbols 

a = radius of curvature 

1 = panel parameter |Eq. (1.17 

b = panel width 

B = panel parameter [Eq. (1.17 

3s = curvature parameter |Eq. (1.17 

a = panel parameter [Eq. (1.17 

D = plate stiffness |EA3/12(1 — v*) 

D = aerodynamic coupling parameter | Eq. (1.17)] 

E = modulus of elasticity 

F = Airy’s stress function 

g, g = structural damping coefficient 

h = thickness 

H = panel rise height (Fig. 1 

k,k = frequency parameter |Eq. (2.4) 

K, Kn = frequency parameter /Eq. (1.14) 

l = axial length 

Kx. 9,4) = aerodynamic loading 

m = longitudinal mode number, number of half 
waves in the axial direction 

M = mach number 

n = circumferential mode number, number of 
half-waves spanwise on the panel, num- 
ber of full-waves around the circumference 
of the cylinder 

V,, Ny, Ney = mid-plane tensile and shear loads 

g = dynamic pressure 

i #. = prestress parameters |Eq. (1.14 

i, = inertia coefficients [Eq. (2.10) 

U,V, 2 = axial, circumferential, and radial displace- 
ments 

u , u = magnitudes of displacements associated with 
mode mn 

U = free stream velocity 
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xy, 2 = Cartesian coordinates 

8 = VM- 1 

4 = aerodynamic damping coefficient |Eq. (2.8) 
6 = solution damping coefficient {Eq. (2.10) 

k = solution frequency parameter |Eq. (2.10 


d = flutter parameter |Eq. (1.14 

v = Poisson’s ratio 

p = density 

? = circumferential coordinate 

x = frequency parameter |Eq. (2.10 

w = frequency (rad/sec 

Q = solution frequency [Eq. (2.10) 
Subscripts 

( ) = air 

( = critical 


material 


Introduction 


eo FLUTTER as a design problem dates back no 
further than the observations of Jordan on the V-2. 
Since that time it has developed to a major design con- 
sideration for any vehicle operating at high dynamic 
pressures by fixing either minimum skin gauge or frame 
The variety 


spacing, over a portion of the structure. 
ratio, 


and importance of secondary effects—aspect 
shape, curvature, edge support, edge loading, pressure 
differential, buckling, etc.—have led to numerous, and 
often contradictory, analytical efforts and experimental 
studies. In spite of these apparent contradictions, 
the need for effective design criteria has led to continued 
consideration of this problem. 

The present investigation was undertaken to deter- 
mine the effects of curvature, a prime factor in any prac- 
tical design. The majority of experimental and theo- 
retical results available relate to flat panels; application 
of these criteria to curved surfaces leads to somewhat 
disappointing requirements. Furthermore, one is in- 
tuitively led to anticipate that the “stiffening due to 
curvature’ should result in considerable increases in 
the natural frequencies, and hence, possibly, in critical 
flutter speeds. Such considerations suggest an ap- 
proach from the standpoint of a correction for the effect 
of slight curvature, using Reissner’s shallow-shell 
theory, on the flutter speed of a flat plate. 


tions of such an approach are reported elsewhere,’ but 


The implica- 


since it does provide an indication of the factors of im- 
portance for the complete cylinder or panels with con- 
siderable curvature, this approach will be initially 
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adopted. Based on these results, a more exact formula- 
tion will be presented, using Gol’denveizer’s equilibrium 
equations for a thin cylindrical shell. 

An alternative approach to the consideration of 
flutter of a cylindrical shell already exists in the litera- 
ture-—-notably the work of Miles,® Leonard and Hedge- 
peth,® and Stepanov.’ With slightly differing structural 
and aerodynamic formulations, these authors consider 
the critical speed in terms of the traveling-wave solution 
for an infinite cylinder without end conditions. This 
paper is not intended to discuss the pros and cons of the 
traveling wave vs. the standing wave as a critical flutter 
solution, nor the differences in traveling-wave solutions 
already presented; hence the title of this paper. Suf- 
fice it to say that there exist cylinders, sufficiently short 
in effective length, that a traveling-wave solution is in- 
admissible on physical grounds. 

Finally, there are already available solutions of the 
standing-wave flutter problem by Stepanov’ and by 
Strack and Holt,* according to an approximate equation 
due to Gol’denveizer.’ Both use differing aerodynamic 
approaches, but are led through approximation to 


SD 
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BS ra 
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Fic. 1. Shallow cylindrical pane’, 


similar differential equations which can be solved in 
exact manner, for arbitrary boundary conditions. The 
approach considered here involves a more exact formu- 
lation of the structural problem and aerodynamics 
which are shown to be similar to either of the above, but 
results in a differential equation which suggests a modal 
solution. Therefore only the ‘‘freely-supported” cylin- 
der will be considered, and the results here should be 
complementary to the above, indicating the limitations 


of each. 


Shallow Shell Formulation 


We begin with the differential equations as derived 
by Reissner! governing finite transverse vibrations of a 
shallow shell described by (see Fig. 1) 


z = 2(x, y) (1.1) 
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In this we assume further that the deflections w(x, y) in 
the direction of the Z-axis are small deviations from a 
constant pre-stressed condition determined by the (ten- 
sile) mid-plane loads N,, V,, N,,. The resulting equa- 
tions, linearized with respect to deflection and mid-plane 
stresses, are given by 
DV4w — Nwr, — Nyy + 2NayWry — 2rrF yy 
—Zyyl'cx + QSryF cy + Puhr = U(x, 9,1) (1.2 
Vik = —Eh (Sp yy + SyWrr — 2WZry'ry) (1.3) 
In the above, F is Airy’s stress function, D the plate 
bending stiffness, A‘ the iterated Laplace operator, p, 
the material density, / the thickness, / the modulus of 
elasticity, and /(x, y, t) the force per unit area acting 
on the shell. 
Following Reissner, we specialize to the cylindrical 
shell with 
2/H = 1 — [y/(6/2)]? (1.4) 
where // is the central rise height, and 0 is the span. 
In terms of the midplane radius, a, 


a? = (b/2)* + (a — H)? (1.5 


or 2aH = (b/2)? (1.6 
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Fic. 3. Frequency vs. longitudinal mode number, m, for example 
cylinder. 
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nd Zrr = Sy = U a 
i rz z ( (1.7) 


s. = —2H/(b/2)? = —1/af 


\We assume the pre-load to be due to a combination of 
pressurization and compression only——i.e., 


Ni = 0 (1.8) 


Further, we take the transverse load as that due to two- 
dimensional, supersonic steady-state flow past the 


evlinder, 


Hx, ¥; t) (—2q/B)w,(x, y, t) (1.9) 


where g is the dynamic pressure and § the reduced 
Mach number. The combination of Eqs. (1.4)—-(1.9) 
reduces Eqs. (1.2) and (1.3) to: 


DAtw — Nw, — N Wy + (1/0) Fr 4+ 


/ JU 


phw;, + (2q B)w () (1.10) 
A‘F (kh/a)w,, (1.11) 


rhis can be compared directly with Hedgepeth’s formu- 
lation for the rectangular panel* (Hedgepeth assumes 
positive compressive mid-plane load). The only dif- 
ference involves the addition of a correction for curva- 
ture involving the stress function. 

Assuming now a solution corresponding to simply- 


supported edges, 
u(x, y, t) W(x) sin (nry/b)e| 


F(x, y, t) 


oy (1.12) 
F(x) sin (nay) /b)e"*'4 


ind reducing to nondimensional terms, 


WY — 2n2W,,"92(1/b)? + niwi(1/b)4W,, + 
RW,” + Rn?x?(1/b)?W, + /a)F,,” + 
x°W, + AW,’ = 0) (1.18) 
FIV — Qn292(1/b)2F,” + 
n'n4(1/b)4F,, = (EhI?/a)W,,” 


where the primes denote differentiation with respect to 
x//) and, following Hedgepeth, 


R, = — N,/° D R, = —N,I° D ’ 


X = 2@i*/@D K? = (phl'/D)wf “4 


We note, in passing, that the solution is decoupled in 
the cross-stream direction; # is only a parameter of 
the solution. At this stage there are several choices 
available, including the exact solution. However, we 
choose the method of modal superposition, relying in 
part on the conservatism of results obtained by Hedge- 
peth for the flat plate. We therefore take, for simply 


Tl 


supported edges fore and aft, 
Sw, sin (marx )| 


W,,(x) 
F(x) = DF» sin (max I) _— 


Substituting into Eqs. (1.10) and (1.11), combining, 
then multiplying by sin (rx //) and integrating, results 
in the set of equations, 


(A r)>> D._w,, = 9 
(1.16) 


m4 — A,m? — By + Cn) mn 


where 
A, = R,/r* — 2n?*(l/b)’ 
B, = (R,/27)n*(l/b)? — n*(l/b)* + K2/ x4 
Can = Cm'*/[(m? + n?(l/b)?|° 17 
’ ) 4 9 . mye 
( = [12(1 — v*)/m4](1/h)*(1/a)’ 
lmr (r*> — m*) r + m odd 
ae 
(0 ry + m even 


We have arrived at a stage where, if a sufficient number 
of m-modes are used, a Galerkin solution approaching 
the exact will result. Reference 3 has shown that for 
thin shells, Reissner’s shallow-shell equation is adequate 
for transverse vibrations of shells corresponding to a 
quarter-cylinder. Further, we shall see that, neglect- 
ing aerodynamic considerations, these equations are 
satisfactory for predicting one of the possible flutter 
modes of the complete cylinder. However, at the 
moment we are principally concerned with determining 
the important parameters in the instability of cylindri- 
cal panels; we therefore pursue a two-mode solution. 

For any value of 7, and values of m and r which differ 
by an odd integer, the two-mode solution will result in 


the following characteristic equation: 


m! — A,m* — B, + C, —(A/r)D, an 
—(A/m4)D,n r§— A,r’ —B,+C 
1.18) 
Noting that Fes —D 1.19) 
and solving for A, we obtain 
A = (44/D mr) [(—m1 + Anm? + By — Cn) X 
(r? — A,r? — B, + C,,)|'" (1.20) 


For fixed values of m, n, and ry the maximum value of \ 


for real values of B,, (and hence w) is then 


A. = [w4(r? — m*)?/Smr| [r? + m? — A, + 
(Cyn — Cmn)/(r? — m?)| (1.21) 


and will occur at a frequency determined by 


B, = (1/2)(m4 + r4 — A,(m? + 7?) + Cran + Cra] 
(1.22) 


Eq. (1.21) corresponds to that derived by Hedgepeth 
for the flat panel, with a correction for curvature. We 
will digress here to consider the magnitude of such a 
correction and the implications with respect to choice 
of modes for the cylindrical panel. For the simply- 
supported rectangular panel, the frequency varies as 
(m?/l? + n?/b*); hence the mode of minimum fre- 
quency corresponds to m 1, m = 1, and for a fixed 
value of m (or m) the frequencies and their differences 
increase uniformly with m (or m). The usual practice 
in aeroelasticity is to choose the lowest order modes 
(both in frequency and mode number) and employ a 
criterion related to frequency ratio squared to fix the 
number of modes required. For the panel flutter 
problem, the modes are uncoupled in m, and Hedge- 
peth’s modal solution follows the stated practice in that 
for the two mode solution he uses m = 1, r = 2 in Eq. 
(1.21), with zero curvature to get (for KR, = 0): 
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Ae = (971/16) [5 + 2n?(1/b)? (1.23) 
For this case, the correction becomes 


Ad, = (324/16)(Con — Cin) 
_ Ont C n*(1/b)?[8 + 5n?(L/b)?} (a0 
16 [4 + nX(/b)2[1 + nl/b42 

To get an order of magnitude comparison, Table I 
presents the results of the two-mode solution for a 
cylindrical panel corresponding to one-quarter of a 
complete cylinder. The dimensions used are those of an 
example to be employed later: length-to-radius ratio, 
l/a = 2; radius-to-thickness ratio, a/h = 500; radius, 
a= 20in.; v = 0.3. The table shows two values of 
\, for the flat panel, the first corresponding to a panel 
with span equal to the width along the circumference, 
the second with span equal to the plan-view width. 
The last column shows the correction for curvature, 
AX, to be applied to the latter value of \,. Values are 
shown for several values of m, and several m and r com- 
binations. In these calculations, C = 6 X 10° was 
used, corresponding to the “effective radius of curva- 
ture’ according to Eq. (1.7). 

The case considered is obviously an extreme example 
in that the ‘‘corrections”’ are several orders of magnitude 
larger than the basic parameter, \,; however, the re- 
sults do illustrate several important points. First, the 
concept of an equivalent flat plate is unreasonable, since 
the minimum flutter speed does not occur in the range 
of minimum mode number; a further explanation of 
this fact follows later. Secondly, even for the two- 
mode analysis it would appear to be a trial and error 
process to locate the modes for minimum flutter speed. 
Finally, while a two-mode solution cannot be considered 
reasonable for an arbitrary pair of coupled modes, the 
results indicate that there may be several modal com- 
binations with nearly the same critical speed, so that 
the particular combination of parameters may be quite 
critical. 

We have thus far considered the problem in purely 
mathematical terms; consider now the physical im- 
plications of the situation. Eq. (1.16) must, for \ = 0, 
yield the equations for determining the natural fre- 
quencies, and, in fact, corresponds to Reissner’s fre- 
quency solution for shallow shells, with appropriate 
notational changes and in slightly different nondimen- 
sional form: 


KR? nn/ tw! = [m? + n?(L/b)?]* + Cm4/ [m? + 
n*(1/b)?|? + (N,1?/9?D)m? + (N,/?/2?D)n*(1/b)? (1.25) 


As has been pointed out by many others, the first 
term on the right represents the effect due to wall 
bending, the second that due to mid-plane stretching, 
and the last two the effects of mid-plane preload. Con- 
sidering the magnitude of the constant, C, in the pre- 
vious example, and neglecting preload, it is apparent 
that stretching effects predominate for low mode num- 
bers. Furthermore, for fixed m, the frequencies will 
always increase with m; whereas, for fixed m, within 
the range of predominant stretching influence, the fre- 
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quencies will decrease with ”; in the range of pre- 
dominant bending effects the frequencies will again in- 
crease with m. Hence a frequency minimum must exist 
for some value of » greater than one for all longitudinal 
mode numbers, m, up to the range where bending ef- 
fects predominate. This leads to the anomalous result 
first pointed out by Arnold and Warburton,‘ that the 
mode of minimum frequency no longer corresponds to 
the mode of minimum mode numbers. 

If we make use of Eq. (1.25) as a definition, we can 


rewrite Eq. (1.16) as 
(Kinn? — K?)t'mn — DX Din = 0 (1.26) 
and the two-mode solution becomes 
As = (Rye — Bg) /2O we (1:27) 
and the frequency at flutter is given by 
K* = (K,,” + Au,*)/2 (1.28) 


Hence we see that the aerodynamic term plays the role 
of a negative stiffness in coupling the zero-airspeed 
modes. Whereas stiffness coupling of a mechanical 
system always drives the frequencies apart, here the 
coupled frequencies are drawn together, until they be- 
come identical. This also indicates the effect on the 
two-mode solution of including additional modes. 
Since coupling occurs only with alternate longitudinal 
modes, the addition of a third mode should stabilize one 
of the two modes without affecting the other directly. 
One would therefore expect a larger value of \,. A 
fourth mode should partially offset the effect of the 
third. Hence fora system with modes well-ordered both 
in mode number and frequency, one could anticipate 
that the two-mode solution is conservative, the three 
mode unconservative, etc. Hedgepeth’s flat panel 
analysis verifies this observation and his four-mode 
analysis approaches, very nearly, the exact. In any 
case, it would appear that an even number of modes 
should be included. 

Next we consider briefly the effects of structural 
damping. In the absence of structural damping, the 
magnitude of the instability would be indicated by the 
imaginary part of the frequency parameter, 


ee => ih. + hans 2 ree — 
(Ren? — Ramada 6 CE28) 


With structural damping included, the potential energy 
terms of Eq. (1.26) would be modified by a factor (1 + 


1Zmn). Presuming equal damping in all modes, Eq. 
(1.29) would become 
K? = (Kyrn? + Kmn?)(1 + tg)/2 + 
t[A2D nr? — (Ken? — Kmn?)?(1 + ig)?/4]'/? (1.30) 
Consider this result in the neighborhood of the previous 
solution, Eq. (1.27), and assume g < 1. Then the 
frequency parameter is given by 
be = Ve a Ce = (Te n* a ea 2" a 2 T 
(Ken? + Kin Qg = (Kin? — Kin*)g!!*]/2 (1.31) 


The flutter frequency itself is apparently little af- 
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TABLE | 


Critical dynamic pressure parameter, Ae, and curva- 
ture correction for a typical quarter-cylindrical panel 
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fected by structural damping. However, at the speed 
which previously predicted the flutter boundary for 
zero damping, the imaginary part of Eq. (1.31) may be 
positive, zero, or negative, depending on the particular 
values of the frequency parameters and damping coef- 
ficient. Should the latter occur, it would obviously be 
a weak instability and would alter the boundary and 
predicted negative damping of the motion at higher 
dynamic pressures very little. It is interesting that 
the possibility of instability due to structural damping 
does exist, a result also noted by Miles’ in his traveling- 
wave solution for the infinite cylinder. However, in 
general, the first term in the imaginary part of Eq. 
(1.31) will predominate, and structural damping, in the 
two-mode solution, will contribute to stability by an 
amount proportional to the sum of the squares of the 
uncoupled modes. This would suggest that in any 
system where the frequencies are not regularly ordered, 
and there are several critical combinations for flutter, 
the modes of lowest frequency will be least effected by 
the presence of structural damping. 

Finally, we mention the effects of aerodynamic damp- 
ing. Borrowing from the discussion of a following sec- 
tion, if we replace the static approximation by a 
quasi-steady form, Eq. (1.9) for the loading is replaced 


by 


No 
~] 


— 2 wy at 
I(x, y, t) = 3 : (x, ee pee w] i 


where y is a coeflicient of order one, dependent on the 

Mach number. Carrying through the analysis, this 

results only in a complex value for the frequency pa- 

rameter expression, in the absence of structural damp- 

ing: 

K? — iK(yd/ UL) (D/ ph)!" = 

(Ken? + Kua*)/2 = (MDa? — (Kyu? — Kusz?*)/4] 
(1.32) 


Again considering the range of the undamped solution, 
Eq. (1.27), 


K? = (K,,? + K,,")/2 + G (1.33) 


where, presuming the frequency is not greatly effected 
by the damping, 


‘ 7 (#4) Ge ~ ==) : ; 
C=, (1.34) 
V 28 \Pnh 2D, 


where p, is the density of air. The first factor is of 
order one, as is D,,,; the second is a relative-density 
parameter, normally quite small. The final factor in- 
volves the product of the sum and difference of the 
squares of the uncoupled frequencies. For a system 
with critical speeds corresponding to high frequency 
modes, this damping contribution may no longer be 
negligible, as is the case for the flat plate. 

To summarize, then, the considerations resulting 
from the simplified analysis of the shallow cylindrical 
panel are as follows. 

1. The concept of an equivalent flat plate is so 
highly conservative as to be unreasonable (Table I). 
The distribution of natural frequencies will not 
be well-ordered with mode number; the minimum will 


occur at some higher value of » [Eq. (1.25)]. 

3. Neglecting structural and aerodynamic damping, 
the dynamic pressure at flutter depends on the differ- 
ence of the squares of the uncoupled frequencies |Eq 
(1.27)|, and the flutter frequency is the root mean 
square of the uncoupled frequencies [Eq. (1.28)]. 

4. In view of the above, it would appear to be ex- 
tremely difficult, if not impossible, to determine before- 
hand the critical range of mode numbers for the flutter 
solution, and, in fact, there may be several ranges, 
nearly equally critical. 

5. Since the frequencies of the cylinder are consider- 
ably higher than those of the flat plate, structural 
and/or aerodynamic damping may be of importance, 
at least in differentiating between flutter modes of 
nearly equal critical speeds [Eqs. (1.31) and (1.33)]. 

It would appear, then, that if a modal approach to 
the solution of the cylindrical flutter problem is to be 
used, the determination of the natural vibration fre- 
quencies should be as accurate as possible, and that the 
formulation should include structural and aerodynamic 
damping. The following section presents a brief deriva- 
tion based on these ideas, followed by a section showing 
the results for a typical case. 


Complete Cylinder Formulation 


The first step in the formulation is to determine 
equations of motion which will adequately describe the 
natural vibration characteristics of the complete 
cylinder. A considerable amount of effort has gone 
into approximate formulations, adequate for particular 
ranges of the number of circumferential modes com- 
pared to thickness-to-radius ratio. In particular. the 
work of Gol’denveizer presents the various equations 
appropriate to the static loading case. The Donnell 
buckling equation is illustrative for large values of the 
circumferential wave number, ». The _ Reissner 
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shell formulation, leading to Eq. (1.25), is equivalent 
to the Donnell equation, modified to include the radial 
In the cylindrical panel flutter solutions 
and §, this equation is further 


inertia terms. 


of references 5, 6, 7, 
modified to include the radial air-forces. 


also considers equations according to Flugge, wherein 


Reference 6 


additional extensional-stress terms are included, and 
references 7 and 8 employ a modified form of the Reiss- 
ner equation, appropriate for shells of medium length 
and nm? < 12 a/h. Thus all of these formulations rely 
on equations of motion which are restrictive on the 


applicable range of m-value, and further, none include 





l—~yp l+yp v l 
Ure + — en Vr -— ~W:— 
Qa2 2a ‘ a 
l+p l—~yp l : 
j ~ tlre =i Ure T Ube — 9 al = V)Urr + 
2a 2 a’ 12a? ( 
h? (2 
12a? 
v, l wT. cd oo 
—U, — —9 (2 — v)Vr2¢ 1Dee w 
2 > “o o.° rh 9 ooo e 
a a? 12a? . a? a’ 
where @ = v/a is the angular coordinate in the circumferential direction, and Y, J, 
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the effects of mid-plane inertia. In view of the appar 
ent importance of exact evaluation of the natural 
frequencies, it would seem that mid-plane inertia ef 
fects could be important, particularly for very low cir 
cumferential Further, it 
desirable that the formulation should not be limited a 
The result must then 


mode numbers. would he 


priori to particular ranges of ». 
be a set of three simultaneous differential equations in 
u 


the three components of displacement, w, radial; 


The appropriate 


9) 


longitudinal; and v, circumferential. 
set, linearized for infinitesimal amplitudes and constant 
mid-plane stresses, and satisfying Gol’denveizer’s com- 


patibility requirements,° is given by: 


— vr [TN. . ) 
a : (v, — a) a 4 0 
Ith a 
| | ; 
; Vode a uw | 
l 7 | v ; 
— VJWrre T W gad 1 V,i Yr) Q (2.1) 
‘ a ; It} 
h? 
Vw — 
12 
l1—yp? |. N, 
| | Nite ie +e) CE ed 
Ih a : 


aud Z are the applied forces 


Replacing the latter by the respective inertia terms, and using the ‘‘freely supported’ solutions, 


ul Um» COS (Marx/1) cos noe’ 
vy = v,, sin (mmx//) sin nde™ (2.2) 
w= W,,, sin (mmx//) cos nde" 
Eqs. (2.1) are reduced to the solution of the characteristic determinant, 
: l—»vp., es 
ee , - 
U = n 2 Mi Vu 
+ (1 — v*)k —fis(1 — v*)un + fll — vu 
(—» ; 
| = 9 gr = i 
ll+y ? OO atecd h? 
| wn = -(2(1 — v)u? + n?] + (2 => ven 
| 2 12a* 12a P 
0 (2.3) 
; ‘ , h* 
— (1 — v?)(A,p° k?) + : 
12a 
) h° 
n = me — gs" H~) 
12a 
h* -~ se ae 9 
— ve - ((2 — v)u'n + 1°] — (1 — v?)(A,u? + t,n?) 
12a? 
+ (1 — v*)n + (1 — p*\R 
wherein p» = mra/l, k®? = pa’w?/E, ni, = N,/Eh, 1, N,/Eh (2.4) 


The solutions of Eq. (2.3) will then yield the charac- 
teristic values k,,,, and hence the natural frequencies 
It should be noted that there are three frequencies 
Considering the eigen- 


WOmn- 
corresponding to each mn pair. 
vectors associated with each frequency, we can classify 
them as predominantly longitudinal, circumferential, or 
radial modes depending on the ratios u,,,/W@ ,, and 
,/Wmn associated with each. Since the linearized un- 


9 
Um 


steady aerodynamic theory provides only a radial com- 
ponent of loading, only the predominantly radial mode 
will be considered, although from subsequent results, 
this may not be sufficient for the low ” modes of large 
l/a shells. The approach here is, however, comparable 
to those employed previously in which mid-plane 
inertias were completely neglected. 
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Further, it should be noted that solutions (2.3) do 
not satisfy boundary conditions which provide for uni- 
form vanishing of the displacements along the node 
lines. In particular, the longitudinal component does 
not vanish at the ends, and the circumferential com- 
ponent does not vanish along longitudinal nodes. 
Hence the present solution is exact only for ring-stiff- 
ened cylinders with massless rings and stringers, if any. 
Eq. (2.1) could be easily modified to include the addi- 
tional inertias; correction for nonzero twisting stiffness 
is not as direct. Here again, however, previous solu- 
tions which ignore all but the radial component in- 
herently contain the same defect. 

We proceed, then, assuming a solution based on the 
predominantly radial mode of a massless ring-stiffened 
cylinder. There remains to be considered the form of 
the aerodynamic loading to be assumed. The features 
of the flow which could be of importance are three- 
dimensional effects, curvature, and unsteadiness; we 
consider these individually. First of all, since for a 
complete cylinder the deflection pattern assumed is 
cyclic and continuous, the aerodynamic loading must 
also be cyclic and continuous. Hence the equations of 
motion will again be uncoupled in , and neglecting 
curvature, the problem is that of a purely supersonic 
or simple planform, for which it has been shown!! that 
two-dimensional air forces provide a very good approxi- 
mation for generalized force calculations. 

As regards curvature, there exists a formal solution, 
due to Randall'* and Holt,'* for the unsteady motion of 
cylindrical shells in a uniform flow, in terms of the in- 
verse Laplace transform of modified Bessel functions 
of the second kind. Miles® and Leonard and Hedge- 
peth® use the steady-state form of this result, appro- 
priate to the traveling-wave problem. Reference 6 
considers the problem with no further approximation, 
but a transcendental equation must be solved which 
tends to obscure any direct effect of curvature. Miles® 
employs a first-order approximation valid in the range 
of his solution, which effectively reduces the aerody- 
namic term to the two-dimensional form. Strack and 
Holt,° in their standing-wave solution, begin with the 
unsteady integral-equation form for the loading, reduce 
to the steady-state form, and employ a linear approxi- 
mation to the resulting kernel function. The loading 
is then given by a two-term approximation, the first 
identical to the two-dimensional form, the second de- 
pendent on radius and circumferential mode number. 
However, this latter term is proportional to displace- 
ment and not to slope, hence the frequency at flutter 
(and the damping ratio) may be effected by curvature, 
but, to this accuracy, the flutter speed is not. Finally, 
it is interesting that Stepanov,’ using the Russian ver- 
sion of two-dimensional piston theory (due to Iliushin), 
derives a critical Mach number for the infinite cylinder, 
which differs by unity from the result of Miles. This 
difference is due solely to the replacement of 6 by .V/ in 
the air force expression. 

Lastly, we wish to include aerodynamic damping. As 
has been mentioned, the formal solution for the un- 
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TABLE 2 
Frequencies of a circular cylinder in radians per sec: 1 = 40 in., 
a = 20in., A = 0.04 in 
m/n 0 a © 2 a 2 
“T [| 9542 5887 3302 1992 1270 #875 y,)6)|US6S.)~*«S 
2 9757 8085 703 4944 3640 2726 ( 367 & 
3 978¢ 9324 8197 6859 6 4551 3054 2554 80 
4 9797 9542 8857 7924 6917 5956 5100 4369 3763 327 
9805 9643 9191 853¢ 775¢ 6952 6179 469 4841 4299 

é 9814 970 38° 8900 8304 7651 IRE 343 74 4 
7 9828 9747 9511 9145 8679 B81 7590 027 648 7 
a 9847 9785 960 8951 ~ 8050 7 8 

? 9875 982¢ 7684 1 880¢ 8413 7998 7 

10 9912 9874 760 9577 33 040 8710 8B 98 7¢ 
ll 9963 9931 78 38 9688 9487 9243 8965 BbE 834 24 
12 | 10028 10002 9926 9802 9634 7430 9195 8938 8665 838 
13} 10112 10090 10027 9923 9783 J611 9413 9194 89¢ 8719 
14] 10216 10198 10145 10058 9941 9797 9629-9443 924 ’ 

1 10343 10328 10284 10212 10113 9992 9850 9693 9524 9347 
16 | 10497 10484 10447 10387 10305 10203 10084 9952 9816 60 
17 | 10679 10669 10638 10588 10519 10435 10336 1022¢ 10107 . 
18 | 10892 10884 10858 10817 10761 10691 10610 10519 10422 31 
19 | 11139 11132 11111 11078 11032 10976 10910 1083¢ 10758 067 
20 | 11421 11415 11399 11372 11336 11291 11239 11181 11119 105 
30 | 16427 16429 16435 16446 16461 16481 16506 1653¢ 1657 661 
40 125398 25403 25416 25439 25471 25511 25561 25 - 7 
m/nf _10 12 l lé 17 18 l 2¢ 

1 643 875 1175 1527 1722 1929 

2 1086 1099 1299 1610 1793 1993 

3 1910 1629 1625 1817 1965 2139 

4 2879 «2362s 2147) «2178s? 2395 

5 3842 3170 2790 2668 2692 2764 

6 4727 3973 ©3481 239 3205 3222 

7 5506 4728 4174 3848 3768 3740 

8 6179 5418 4840 4463 4351 4290 

9 6758 6039 467 5067 4937 4852 

10] 7262 6599 6051 5650 511 414 

11] 7706 7106 6596 6209 607¢ 5969 

12] 8105 7570 7105 6745 6612 6512 

13] 8475 8004 7588 7260 138 7045 

14 8826 841¢ 805¢ 7760 é 756° 
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1¢ 9508 9210 8944 8734 8603 
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18 | 10216 10014 9838 9708 7643 

19 | 10593 10433 10299 10207 174 

20} 10991 10871 10775 10719 10718 

30 | 16663 16783 16935 17124 17354 

40 | 25855 26060 26305 26592 692 








steady loading on a cylindrical shell, with cyclic defor- 
mation both in time and space, is available, but not 
in particularly usable form. Alternatively, the above 
discussion implies that the two-dimensional form should 
provide a good approximation for the in-phase com- 
ponent. This would suggest use of two-dimensional 
unsteady theory as a first approximation to including 
unsteady effects. We will go one step further here and 
make use of a form correct only to first order in fre- 
quency. At high supersonic speeds, this is equivalent 
to piston theory,'* which is the exact result. 


Therefore, the expression to be used at low supersonic 


speeds (and WV > 72) is” 


K(x, y, t) = 


2q 3 M? -—2 Wy ] Wi\ r=0 


a 


8 M?—1l M@-—-i1l 

while for sufficiently high supersonic speeds, the expres- 

sion becomes 
I(x, y, t) — (2¢g/.M)[w, + (1/U)u, | 2.6) 


Making use of the fact that the leading edge of the 
shell is restrained and the time variation is simple- 
harmonic, these can be combined into 


—2 wy ; e 
U(x, y, t) = t(w, + i— w)e™ (2.7) 
‘ 8 


where g/8 — g/M for M > MM, and 
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M?-2 V2<M< M) 
y= (MP1 


(2.8) 
] l M> us 


where the value of .\/; is left to the discretion of the 
user. 

It remains, then, to derive the equations of motion, 
employing the predominantly radial modes resulting 
from Egs. (2.1)-(2.4) and the loading from Eq. (2.7). 
Either from Lagrange’s equation, or by analogy with 
the plate equations, these become 


jf QR + a iia X7 len 


2qa* 
DarWm = 0 (2.9 
BElh “ 


1 + (tmn/Wmn)? + (Omn/Wmn)? (2.10) 


where 7’, 


; 2ga? 1/2 ¥ (=) 2 ; 

= Kk? — «(- (2.01 
_— a V8\ pmh 
k = Qa(p/F)'!2(1 + 76) (2.12) 


( and 6 are the solution frequency and damping ratio, 
and the remaining symbols have been defined pre- 
viously. The 7’, terms are the corrections to the 
generalized mass which come about due to coupling of 
the in-plane modes with the predominantly radial 
mode. 

There will now be a set of solution frequencies equal 
in number to the number of modes assumed, each with 
its associated damping ratio. The stability boundary 
will be fixed by the vanishing of the damping ratio for 


one of the modes. For zero airspeed, we will have 


Kmn = Rinn(I =e 12 mn/ 2) (2.13) 


Hence a comparison of the solution damping ratio with 
one-half the assumed structural damping will provide 
an indication of the change in stability of any mode due 
to the influence of the external flow. 


Sample Calculation 


The following calculation was undertaken to deter- 
mine the modal ranges of importance in the standing- 
wave flutter solution. In that the structural represen- 
tation is unlimited with regard to mode number, and 
the aerodynamics employed appear to be a reasonable 
approximation to first order in curvature and frequency, 
the present approach should provide an improved solu- 
tion to the problem, limited only by the adequacy of 
the number of modes considered and the assumption of 
a single, primarily radial mode for each mode number 
combination. To the extent that these calculations 
may be considered typical, they should provide an indi- 
cation of the need for including the various refinements 
considered, the number of modes required in a modal 
solution, and the adequacy of the criteria for modal 
selection based on the two-mode solution. 

The particular values of the parameters chosen are 
aluminum cylinder of 20-in. radius, 40-in. long, with a 
thickness of four-hundredths of aninch. The following 
numerical values were employed in the analysis: 
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kb 10° Ib in. p 0.260 X 
10—* Ib-sec?/in.4 
Nr = Ny 0 v = 0.3 (3.1) 
l/ja = 2.0 h/a 2X 10 
a = 20)in. 


Based on these data, the natural frequencies of the 
cylinder (in rad/sec.) were computed for mode number 
Table 2 presents 


combinations up to m = 40, » = 40 
a reduced tabulation of these results, which are also 
shown graphically in Figs. 3 and 4. 

While Table 2 contains a wealth of data, several 
points are quite apparent. The mode of minimum fre 
quency, for example, occurs for nm = 8, m 1, while for 
n = O the lowest frequency is an order of magnitude 
higher. Further, the first 20 modes for n = 0 are con 
tained within a frequency band of 2,000 rad, sec., while 
for successively higher modes the spacing increases 
This is particularly significant in the flutter solution, 
wherein the modes are uncoupled spanwise 1.e., we 
consider the solution for each value of n. According 
to the two-degree-of freedom results, the dynamic pres 
sure at flutter will be directly related to difference in 
squares of the frequencies; the proximity of the neigh 
boring modes provides an indication of the number of 
modes which will couple strongly 

Regarding the difference of the squares of the fre 
quencies as the product of the sum and difference, the 
susceptibility to fiutter is partly indicated by the slopes 
Ow/Om in Fig. 3, or by differencing in Table 2. It is 
apparent that there exist two distinct ranges of mini- 
mum frequency difference, the first for low at an 
intermediate value of m. As increases, the minimum 
difference and the value of m at which it occurs both 
increase, the former at a considerable rate. At a suf- 
ficiently large value of 1, this minimum difference no 
longer occurs, and Ow/Om increases monotonically. 
These minima all occur at about the same frequency, 
in this case in the neighborhood of 8,000 to 9,500 rad / 
sec. The second range of minima occurs for large 
However, in this case the 
Hence, we 


values of m and small m. 
average frequency continuously increases. 
may anticipate a minimum in difference of the squares 


of the frequency for some finite value of ”. These 


W X 10-3 (RAD/SEC) 





° nm 
10 15 20 








Frequency vs. circumferential mode number, 7, for 


Fic. 4. 
example cylinder. 
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TABLE 3 


Parameters for estimation of critical flutter modes 
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effects are summarized in Table 3, which presents, for 
several values of m, the minimum frequency difference 
®, — Wm, the longitudinal mode numbers m, r for which 
the minimum occurs, the average frequency at which it 
occurs, (w, + w,,)/2, and the actual difference in the 
squares of the frequency parameters, k,” — k,,” |see 
Eq. (2.4)|. These are given for both ranges men- 
tioned above, the first denoted low n, the second high 
n. While the table is based on minimum frequency 
difference rather than minimum difference of the 
squares, the results would be very nearly identical 
According to this evaluation, there are two sets ot 
mode numbers for which the difference in k® is very 
nearly the same: » = 0,m = 4,r = SJandn 16-20, 
m = 1,r = 2; however, they are considerably different 
in nature. In the first case, since the minimum occurs 
at an intermediate value of m and a considerable num- 
ber of modes are closely spaced in the x = 0 frequency 
spectrum, it may be anticipated that a large number of 
modes will be required in the flutter solution to obtain 
convergence. For the second case, the minimum occurs 
for the minimum m and the higher modes are well- 
spaced; hence relatively few modes should be re- 
quired. Furthermore, using the average frequency as 
an approximation to the root mean square frequency, 
the first case should be considerably more influenced by 
damping, either structural or aerodynamic, than the 


second. 

The flutter computations were carried out in accord- 
ance with these considerations. For low ” (n = 0, 1), 
the initial evaluation employed m = 1-12. After 


locating the flutter range, 16- and 20-mode calculations 
were employed to verify the convergence of the modal 
approach. For the high ” range, four modes were used 
to determine the flutter range, with six- and eight- 
mode calculations to check the accuracy. Two values 
of structural damping were assumed, equal for all 





modes, g 0.001 and g = 0.01. The results are pre- 
sented in Figs. 5 to 7, where the damping ratio, 6, is 
plotted vs. dynamic pressure divided by reduced Mach 
number, g/@, for the mode which first exhibits instabil- 
ity. In all cases, the flutter frequency is very nearly 
equal to the root mean square of the frequencies of the 
modes which couple directly to produce the instability 

Figs. 5 and 6 present the damping of the critical 
modes for the basie calculations at g 0.001. The 
low n case was initially performed using 12 modes for 
several values of nm. The result shown in Fig. 5a corre- 
sponds to m = 0; calculations at considerably higher 
values of g/8 did not show an instability for larger 
values of nm. Asa check on the 12-mode solution, a 16- 
mode calculation was run. This resulted in an increase 
of about 5.5% in flutter g/8 as shown in Fig. 3b; asa 
point of interest, the mode which exhibited the insta- 
bility was shifted from the fourth to the fifth. A 20- 
mode calculation verified the 16-mode results 

In Fig. 6a are given the results of a four-mode calcu- 
lation for the high m case; only three of the modes 
Actually, 
modes of higher and lower » become unstable at dy- 


which first exhibit instability are presented 


namic pressures just slightly higher than those shown 
A six-mode calculation led to a 4.6% increase in dy- 
namic pressure at flutter, as shown in Fig. 6b, and the 
result was confirmed by an eight-mode solution. In all 
cases the m = | mode was critical. 

The present results indicate that the initial instability 
for very low structural damping occurs in the mode » 
Qatg/8 = 5.4; the next occurs at g/8 15.1 for » 
16. However, the instability at 0 is relatively 
weak by comparison (note the scale difference between 
Figs. 5 and 6) and occurs at a considerably higher fre- 
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Fic. 6. Damping ratio vs. reduced dynamic pressure for 
high n 


mode would be considerably more affected by structural 
damping. 

Fig. 7 presents the results for the two cases for g = 
0.01. A 16-mode analysis was used for m = 0 and the 
seventh mode proved critical; for the high ” case four 
modes were used. It is apparent that the critical con- 
dition for nm = 0 has been nearly tripled, while the high 
n condition is but little affected; both are now of the 
same magnitude, although the high 7 results should be 
corrected with a six-mode solution. 

The effect of aerodynamic damping was not included 
in the previous calculations. However, for the com- 
puted damping ratio magnitudes (6 < 1), and consider- 
ing the magnitude of the aerodynamic damping coef- 
ficient (approximately 10~* for the critical speeds at 
sea level), Eq. (2.1) can be closely approximated by 


iy [2qa* \}!? ( Pal ) . Pad ’ — 
(1 + #8) (3.2 
dees (=) —— .—* ) GH 


Since | x| = |, the effect of aerodynamic damping can 
be estimated by adding one-half the aerodynamic damp- 
ing coefficient to the computed damping ratio. From a 
consideration of Figs. 5-7, it would appear that aero- 
dynamic damping could delay the low x instability to 
In the case of 


II? 


somewhat higher dynamic pressures. 
the high m modes, aerodynamic damping should have 
relatively little effect in shifting the stability bound- 
aries. 

One further comment should be made concerning the 
accuracy of these calculations. The program was 
based on a value of 7°,,,, equal to unity [see Eq. (2.10) ]. 
Now, if Eq. (2.9) is divided by T,,,,, the reciprocal may 
be considered an aerodynamic effectiveness factor. In 
particular (see reference 16), for the range m > 2, T,, = 
(1 + n~*), therefore, the computed dynamic pressures 
could be corrected by multiplying by the reciprocal, 
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Tm, '. Obviously, for the high » flutter cases, this 
factor is negligible. For the low n cases, the following 


can be determined from reference 16: 


Tyo = 1.087 1x9 1.011 
Ty, = 1.534 ZT» = 1.0604 
Ty. = 1.247 Tx. = 1.090 
Ti3 = 1.128 1s; = 1.079 
74 = 1.070 Tx = 1.058 
Ti35 = 1.044 To = 1.042 


Considering the results obtained, the effect of the cor 
rection should merely reduce the coupling of the lowest 
one or two m modes for low n. Since these modes are 
reasonably remote from the critical m mode, the results 
presented should not be greatly affected. Again, a 
check calculation would be desirable. 


Discussion and Conclusions 


The procedure presented is admittedly painstaking 
and time-consuming. Furthermore, the nature of thc 
solution is more nearly applicable to check stability of 
a given configuration than to determination of required 
parameters for preliminary design. Finally, the prob 
lem considered is idealized in that the mass of frames 
and stringers has been neglected. However, this last 
consideration is no more restrictive than previously 
mentioned solutions, the aerodynamics assumed are 
considered equivalent, and the structural equations in 
clude both the exact coupling forms for all ” and the 
effects of mid-plane inertia. The results should there 
fore provide an evaluation of the range of applicability 
of other solutions. <A singular advantage of the non- 
modal solution lies in the applicability to other than 
simply-supported boundary conditions. 

Reviewing the results obtained, there are two ranges 
of critical modes for the cylinder. The first corre- 
sponds to low values of the circumferential mode number 
n, but relatively high values of the longitudinal mode 
number, m; this will be referred to as membrane-type 
flutter. 
values of m, and in the latter respect is similar to plate- 
type instability. While it would be desirable to discuss 


The second occurs for high values of 7 and low 
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Fic. 7. Damping ratio vs. reduced dynamic pressure for moderate 
structural damping. 
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the problem in terms of the flutter boundary as a func- 
tion of n, the computations involved for nearly exact 
solutions for all m appeared prohibitive, and the follow- 
ing are conclusions concerning the two types of stability 
based only on the calculations presented. 

The membrane-type instability is critical for n = 0; 
this is indicated by the two-degree-of-freedom criterion 
refer to Table 3) and also by the calculations per- 
formed for low 7. The instability occurs for n 0 in 
the range of m predictable on the basis of minimum 
frequency difference. For low structural damping, and 
the particular radius, length, and thickness, this is the 
critical instability of the evlinder. However, the criti- 
cal value of dynamic pressure is severely effected by 
structural damping, and aerodynamic damping may be 
of importance. The frequency spectrum for n Q is 
very closely spaced and a modal solution will require a 
large number of longitudinal modes to obtain an accu- 
rate evaluation; the alternative of an exact solution of 
the differential equations would appear desirable with 
the condition that mid-plane inertia be included. 

That this latter condition is essential is reasoned as 
follows. Neglect of the mid-plane inertias for a thin 
shell leads to the Donnell equation with radial inertia 
loads.* This is exactly equivalent to the shallow-shell 
vibration theory derived by Reissner,' which for n = 0 
reduces to |see Eqs. (1.25), (1.14), and (2.4) | 


(ma)! hk\? fa\? 
knw = 1 + (4.1) 
12(1 — v?) (") (") 


Had this expression been used, it would have resulted 
in a minimum difference in the square of the frequency 
parameters for longitudinal mode numbers | and 2, 
with a value kw? — Ry" 2.43 X& 10°-° as compared 
with the minimum value kj” — kw? = 1.67 & 10 
actually calculated (see Table 3). While the combined 
effect of a change in magnitude of the difference and in 
mode number at which the minimum flutter speed occurs 
is difficult to assess, we may compare the results for 
the present case with that given in reference 8, which 
infers, for n = 0, freely supported end conditions and 


no structural damping, 
q/B = |l5E/( — v*)|(h/))* (4.2) 


For the values of the present case this predicts g/8 = 
0.165. The sole differences between this result and 
that calculated herein (¢/6 = 5.4) lie in a relatively 
minor difference in structural coupling, and the inclu- 
sion in the present calculation of mid-plane inertias and 
structural damping, g = 0.001. To eliminate the latter 
difference, a zero damping case was run for several 
values of the dynamic pressure, and these indicated a 
value g/8 = 4.5. 

For the panel-type flutter,* the instabilities occur at 
sufficiently high circumferential mode numbers () that 


* Following preparation of this paper, the work of Shulman 
was brought to the author's attention. Reference 17 treats the 
high » case only, also by a modal approach, but without damping 
Results and conclusions are similar to those presented here, within 


these limitations. 


the frequencies are well ordered, m = | is the critical 
mode, and the critical speed can be well predicted with 
a modal analysis using relatively few modes. For this 
range of m, the frequencies are given with good accuracy 


by shallow-shell theory [Eq. (1.25)], and the range of 
critical 7 can be predicted on the basis of the two-mode 


solution 
0/On(k2,”? — k,7) = O }.3) 


where for the complete cylinder we replace 6 with 7a. 
The result is, for » > mza/1/, 


n. = (6001 — v*)r?]6(a/1)"8(a/h) 4.4) 


where 7 must be an integer. For the present case this 
yields n, = 18, which is in reasonable agreement with 
the calculations given here. The undamped two-mode 
solution of Eq. (2.9) for n 18 gives g/8 = 9.69, a 
conservative estimate of the value of the calculated 
flutter speed. 

Eq. (4.4) implies a functional dependence 


qg/B ™~ E(h/a)?/*(h/1)> 1.5) 


subject to the requirement of integer 

For low structural damping and the particular case 
considered, the panel-type instability is not critical. 
However, the dynamic pressure at flutter is relatively 
unaffected by structural damping, and for a value of 
damping g = 0.01, the panel-type and membrane-type 
instabilities are of nearly equal importance for the case 
presented. Further, the negative damping ratio for 
panel-type flutter greatly exceeds that for m = 0 at 
speeds above critical. 

It is to be noted that the present calculation resulted 
in a finite value of critical speed at a specific m rather 
than varying as n~‘ as predicted in references 7 and 8 
using the Gol’denveizer equation for ‘“‘medium’”’ shells. 
This is not surprising, since the inequality n> < 12 
(a/h) = 1,730 is not well satisfied for the range of the 
critical mode. Reference 8 suggests that the flutter 
may occur in the minimum energy mode, which is also 
the mode of minimum frequency. From Eq. (1.25) 
this may be determined as 


ny & [121 — v?)r*]'8(a/1)'/*(a/h)'3 (4.5) 


which here gives m) = 8 in agreement with the calcula- 
tions, whereas the flutter mode was n = 16. 

Throughout the preceding discussion and analysis, 
the two-mode solution was used as a guide. In no in- 
stance were the indications of the two-mode solution 
contradicted as regards the importance of frequency 
differences (and therefore mid-plane inertia effects), 
the selection of critical mode ranges, the effect of addi- 
tional modes, or the effect of structural damping. 

It would appear that the modal approach outlined is 
entirely satisfactory for the determination of the panel- 
type flutter which occurs at high circumferential mode 
numbers. In fact, the frequency computation for these 
modes according to Eq. (2.1) could be replaced by Eq. 
(1.25) with no loss in accuracy. A conservative esti- 
mate for this type of instability may be gained from the 
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two-mode solution for m = 1, 2 and n for several values 
near m, as determined above. The addition of higher m 
modes in pairs results in a more accurate and always 
conservative result. Analysis with an odd number of 
modes will yield an unconservative result. However, 
as has been seen, for very low structural damping this 
may not be the critical instability of the cylinder. 

For the membrane-type flutter, the modal approach 
is not entirely satisfactory due to the large number of 
modes required to obtain convergence. This difficulty 
would be further emphasized for smaller thickness-to- 
radius ratios. While an exact solution would be prefer- 
able, mid-plane inertias must be included. Conver- 
ence of the modal solution can be obtained with a suf- 
ficient number of modes, but conservatism of any 
simplified calculation cannot be guaranteed. This 
brings up one further difficulty, common to all the solu- 
tions mentioned. In treating the problem in terms of a 
single mode or equation for each value of 7, we have as- 
sumed a single predominantly radial mode. However, 
for low values of m and length-to-radius ratios slightly 
greater than those considered here, it becomes difficult 
to distinguish the predominantly radial mode, and if 
one exists it may not be the mode of minimum fre- 
quency. For example, for » = 0 the circumferential 
displacement is decoupled from the longitudinal and 
radial motions and may be neglected. The two re- 
maining modes, their displacement ratios and _ fre- 
quency ratios are then given as follows (see reference 
16): 


l/a (Wy Uo)? (Wyo /Uw) 2) 0°” /an0' 1) 
2.0 —3.406 0.2936 1.74 
3.0 —1.165 0.8580 1.89 
3.5 —0.7025 1.422 1.39 
4.0 —0.4760 2.101 1.48 
6.0 —0.2073 4.825 2.07 
10.0 —0.10385 9.661 3.37 


It would appear that an accurate analysis for cylin- 
ders with //a > 3 would require consideration of both 
modes in a modal analysis, or equivalently, both 
equations of motion in an exact solution. For n = 1, 
the situation is further complicated by the coupling of 
circumferential motion and the apparent need to in- 
clude all three modes of motion in the analysis. For 
n > 2, the mode of minimum frequency always corre- 
sponds to the mode of maximum radial displacement, 
and the frequency ratios of the higher modes would ap- 
pear to preclude their importance. 

Finally, it should be mentioned that the form of 
structural damping, equal hysteresis damping in all 
modes, has not been substantiated by any tests. In the 
experimental results given in reference 3, it was ob- 
served that the damping was more nearly of the viscous 
variety, with gw = 12 m rad/sec. for all modes. In 
view of the demonstrated importance of damping, 
further experimental evidence would be very desirable. 
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Appendix 


Throughout the above analysis, a form of structural 
dissipation was assumed which is proportional to the 
strain energy, but in phase with velocity. This repre- 
sentation is only an approximation, but has been found 
to be satisfactory for nearly sinusoidal motion in un- 
coupled modes, and is commonly employed in flutter 
analysis. 

In the present application, this imposes certain 
limitations on the interpretation of the numerical re- 
sults, principally due to the form of the assumed solu- 
tion as a frequency: 


w = A1 + 78)\ 


int _  ,— 8M im (A.1) 
= ¢ é 


(Continued on page 961) 
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Leading-Edge Separation of Laminar 


Boundary Layers in Supersonic Flow 


Tt 


W. B. BROWER, JR.* 
Rensselaer Polytechnic Institule 


Summary 


A brief description of the flow field is given for the interaction 
of shock wave and laminar boundary layer on a compression 
\ special sub-case—that of leading- 
is analyzed by extension of Chapman's 
laminar mixing-layer theory Results are tabulated for ranges 
of Mach number 1 < J, < 10, and compression-corner angle 
0<8< 45°. A limited region of possible leading-edge laminar 
separation with an attached leading-edge shock (or in certain 
cases an expansion) followed by a second shock due to the re- 
Comparison with existing 


corner in supersonic flow 
edge laminar separation 


attachment flow is found to exist. 
experimental data is found to be satisfactory in several cases. 


Symbols 
VW Mach number 
p = static pressure 
u = flow speed 
tis = Chapman’s speed ratio for laminar mixing 
B = compression-corner angle 
1 = ratio of specific heats 
6 = flow deflection across shock 
4 = shock inclination with respect to velocity vector ahead 
of shock 
g static-pressure ratio across a shock wave 
Subscripts 
1 = denotes shock 
0 = free-stream conditions 
1 = conditions behind leading-edge shock 
2 = conditions behind reattachment shock 
. = conditions on separation streamline 
min = minimum value 
max = maximum value 
Superscript 
: = conditions for reattachment shock at limiting case of 


weak shock 


Introduction 


See PROBLEM OF INTERACTION between shock waves 
and boundary layers has recently received much 
attention—e.g., by Chapman, Kuehn, and Larsen,' and 
by Gadd, Holder, and Regan.’ It is clear in general 
that the interaction is a very complicated phenomenon 
in which the boundary layer may separate or, if orig- 
inally laminar, become turbulent. The status of the 
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theory does not yet permit complete analytical pre- 
diction of the conditions under which separation or 
transition will occur. There is, however, one special 
that of leading-edge separation of a laminar 
which 


case 
boundary layer from a compression corner 
appears to be amenable to a rather simple analysis. 
The basic flow configuration is shown in Fig. 1. If 
leading-edge laminar separation obtains, a mixing- 
layer joins the main stream to a pocket of recirculating 
gas lying between the separation streamline and the 
The ensuing flow deviations create a double 
shock system. The leading-edge shock is caused by 
the deviation of the flow due to the mixing-layer. 
Downstream of this shock, not too near the 
corner, the streamlines are straight and parallel to the 
outer edge of the mixing layer. As the separation 
streamline curves into the reattachment point a system 


body. 


and 


of compression waves is generated; these waves coalesce 
and eventually form a reattachment shock. 
Experimental results, presented by Chapman, Kuehn, 
and Larsen in their Fig. 27, demonstrate that the static 
pressure rise ~2/p, across the reattachment shock is 
independent of Reynolds number and depends only on 
the inviscid-flow Mach after reattach- 
ment. In the same paper they present an analysis, 
based on Chapman’s laminar mixing-layer theory,’ 
which gives the following result for the pressure ratio. 


number Ms, 


= l T (Y — 1) M2 2(1 —_ Ux") ¥/(y¥- 1) (1) 


7 1+ (y — 1)M,*/2 


Str 


It is perhaps desirable to emphasize that this relation is 
entirely independent of the shock relation for the 
pressure ratio £». 

Their analysis is based on the following assumptions: 

1. The static pressure along the straight part of the 
separation streamline is equal to the inviscid-flow 
pressure /y. 

2. Along the curved part of the separation stream- 
line the gas is compressed isentropically to reattachment, 
at which point the flow is stagnated at pressure py». 
Everywhere on the body downstream of reattachment 
and outward from the boundary layer into the inviscid 
flow behind the reattachment shock, the static pressure is 
also p2, since the pressure gradient normal to the sur- 
face is presumed negligible. 

3. In the mixing-layer, similar velocity profiles 
exist; therefore from the Chapman mixing theory the 
ratio of the separation streamline flow-speed ux to 
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Fic. 1. Leading-edge laminar separation 


the inviscid flow speed i; is the constant valuet 


fie = ta /t — 0.987 


1. Perfect gas relations apply. 

5. Heat transfer at the wall is negligible, and the 
Prandtl number is unity. Thus the Busemann energy 
integral can be used to relate the mixing-layer tempera- 
ture to the inviscid-flow total temperature. 

In a series of experiments on a variety of shapes with 
leading-edge separation, Chapman ef a/. (see their Figs. 
27 and 28) have confirmed the validity of Eq. (1). For 
emphasis we note again that Reynolds number plays 
no role in this flow. This of course is not to say that 
effects of viscosity are negligible. 


Extended Theory 


Since for leading-edge separation the reattachment 
pressure is not a function of Reynolds number, it is 
reasonable to expect that inviscid-flow relations, along 
with Eq. (1), should be capable of completely deter- 
mining the flow field exterior to the separation region. 
Referring to Fig. 1, it is assumed that the mixing-layer 
causes the inviscid flow to deflect through an angle 
6, across the leading-edge shock; that is, 6; defines the 
exterior of the laminar-mixing region. Then across 
the reattachment shock there must be an additional 
deflection such that their sum equals the inclination 
of the downstream wall. Across the two shocks the 
usual gas dynamics relations are presumed to hold. 
For the present purpose the appropriate equations are 


- | (2) 


(y — 1). sin? 044. + 2 |!" 
(y -— 1) 


(y + IM? 


cot 6,4; = tan ssl 5, V2 sin? 6 1) 
_ - t - if i+1 -— 


l 
M4 _ 
sin (0;41 


x 
6:41) 


: (3) 
27M? sin? 6:41 — 
where 7 = 0, | denote the leading-edge and reattach- 
ment shocks respectively, and 
t For separation originating at a point of zero boundary-layer 
thickness, Chapman’ has shown that is is essentially inde- 
pendent of the choice of viscosity-temperature function. 
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g = [2y\J/,* sin? — (y — 1)\/(y + 1) (4 


. 
In addition we have the kinematical condition 


6» = p— 6; (2) 


Eqs. (1) through (5) are a set of seven simultaneous 


equations in eight variables, with 


parameters y, Wo, 8. The possibility of multivalued 
solutions is obviated by the additional stipulation that 


transcendental 


for a given flow deflection, the solution corresponding 
to the smallest shock angle is to be preferred. This is 
consistent with experimental observations of weak shock 


Waves. 


Method of Solution 


Due to the transcendental nature of the equations 
explicit forms of the solution are not possible. Use is 
made of a double-iterative numerical solution. The 
actual computations for any extensive range of param- 
eters can be handled only by high-speed digital com 
puter. The results presented herein were obtained on 
AVCO's IBM 704 computer. 


Results of Computations 


It is of interest that in the numerical solution of the 
equations it became apparent that not every combina- 
tion of .\/) and @ yields a solution. <A full discussion 
of the analysis to determine the limiting cases is pre- 
sented in the Appendix. There it is shown that for 
the range of .1/) considered (1 < J) < 10), there is a 
region with an upper limit and nonzero lower limit on 
3 which can sustain leading-edge separation. This is 
plotted in Fig. 2. According to the present theory it is 
remarkable that for Wo < 1.27, leading-edge separation 
is not possible for any value of 8. It is possible that the 
limited range for which solutions exist at l/) = 1.5. 
accounts for the experimental difficulties alluded to by 
Gadd et al. in a footnote, p. 239 of reference 2. 

The numerical resultst for those parameters calcu- 
lated are presented in their entirety in Table 1 at 
the end of this paper. Plots are given in Figs. 3-5 


+ Acknowledgment is made to AVCO staff members Dr. Mur- 
ray Mannos for advice in the analysis of the numerical solution 
and Mr. Leo Morissette who prepared the program and carried 


out the machine computations. 
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LEADING-EDGE 


4 the separation region angle 6,, the leading-edge 
shock angle 6;, and the pressure-ratios p/ po and p2/ po 
jor the two shocks. These four quantities are the 
simplest to measure in a wind-tunnel experiment (see 


references | and 5 for a discussion of techniques). 


Comparison with Experimentt 

Chapman! presents sufficient data on two cases of 
leading-edge separation to permit comparison with the 
present calculations. For the first, designated herein 
as Case A,** data are taken from Figs. 26c and 28, and 














TABLE 1 
ase A ase f 
M - MM - 7 
os 25 (3 = 25 
——— S ee 
eor i Theor Exp 't 
) fa o | 
i 8 4.0 
8 ; 
i 4+U, 5 34 
My L~2o 3 - = 
p,/p, 59 »4 + -? 
< i 
p,/p = - 2.4 2.35 
l 
95/p -- -- +.50 37 


jor Case B from Fig. 4a. 

The limit of experimental accuracy is not known. 
There may also be errors introduced by measuring 
angles directly from the photographs as reproduced 
and published. In any case, in these two examples 
(see Table 1) the comparison is satisfactory, considering 
the simplicity of the theory. 


Discussion 

The simple theory of this paper allows prediction of 
the exterior flow conditions for leading-edge laminar 
separation on a compression corner at supersonic 
speed. It delineates a limited region for parameters 
VU and 8 outside of which leading-edge separation does 
not appear to be possible. 

A nice verification of a lower limit on the compres- 
sion corner angle to sustain leading-edge separation is 
given in Fig. 30b of reference 1, in which it is found 
that for 8 = 10° at My = 2, separation takes place 
downstream of the leading-edge; its position and the 
pressure distribution depend on Reynolds number. As 
Reynolds number increases, the separation point moves 
further downstream and the pressure rise at separation 
decreases. 

Experimental verification of 8,,., is not known for pla- 
nar flow, but in axisymmetric flow over a flanged cylin- 


t The author would appreciate receipt of experimental results 
bearing on this problem. 

** Case A is actually a 20° compression corner inclined at 
—5° to the stream. See the Discussion for conditions under 
which the analysis is applicable. 
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Fic. 3. Flow deflection across leading-edge shock 


der, Lee® finds that above a certain value of 6 the flow 
character changes radically, and depends upon the height 
of the compression-corner base (which is necessarily 
finite for wind-tunnel test purposes) rather than the 
angle 8. His results cannot be compared directly with 
this present work which is restricted to planar flow 
It appears, however, that in this flow regime detached 
shocks appear. Caution must be observed in use of 
Fig. 2, since it is not certain that a case lying in the 
region between 6,,;,, and 8,,,, will automatically admit 
of leading-edge separation. In Fig. 2la of reference |, 
for a 20° compression corner at J, = 
occurs slightly downstream of the leading-edge and 
This seeming anom- 


3, separation 


is Reynolds number-dependent. 
aly can be clarified only by additional testing 


Extension to Nonhorizontal Upstream Walls 

The computed results may be used directly to include 
compression corners for which the upstream wall is 
inclined to the free-stream direction. The only re- 
quirement is that its inclination to the free-stream be 
sufficiently less than the 6, corresponding to 8 (always 
measured with respect to free stream) so as not to dis- 
turb the effective condition—implicit in Chapman's 
mixing-layer theory—that the surface be infinitely far 
from the separation streamline. 


Appendix 


Limiting Cases of Leading-Edge Separation 

It is known from experiment that not every combina- 
tion of 8B and Mp, will yield a leading-edge separated 
flow. The question is raised if the analysis based on 
Chapman’s simple theory predicts limiting cases 
It is shown that the theory yields, for given Jo, upper 
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TABLE 2 
ee 7 Leading-Edge Laminar Separation in Supers« mic Flow = 

M 

o | A] & 5 M | Pi/Po 8, g, Ma | Po/Po 

1.5 | 10°] 45.48°] 3.17! 1.39 |] 1.17 | 57.0%] 6.87] 1.12] 1.66 

2.0 115 | 34.21 4.90] 1.82 1-31 | 43.50 1.10 | 1.47 2.19 

20 | 40.80 | 11.33] 1.59 1.83 | 49.40 8.67 | 1.28 2.80 

2.5 115 | 25.99 3-22 1 2.36 1.23 | 35.24] 11.78 | 1.90 2.47 

20 | 30.66 Beta | 2.3% Sa7s § Seekk | 23688 | 1.72 3.24 

25 | 36-51 | 14.60 | 1.89 2.41 | 42.17] 10.40 | 1.52 4.12 

30 | 44.83 | 21.45] 1.58 3.46 | 49.88 8.55 | 1.27 5.29 

3.0 | 15 | 21.43 2.77.1 2.86 1.24 | 30.55] 12.23 | 2.28 2.85 

20 | 25.54 7.92] 2.61 1.79 | 32.72 | 12.08 | 2.09 3.85 

25 | 30-48 | 13.26 | 2.34 2.5% | 35.49 | 11.74 | 1.88 5205 

30 | 36.59 | 18.99 | 2.05 3-56 | 39.48 | 11.01 | 1.65 644 

35 | 45.20 | 25.76 | 1.67 5.12 | 47.05 9.24 | 1.35 8.09 

3-5 715 | 18.56 2.84 | 3.33 1.28 | 27.31 | 12.16 | 2.65 3.26 

20 | 22.43 7077 | 320% 1-91 1 29.22 | 12.23 4 2.42 4.58 

25 | 27.02 | 12.82 | 2.7% 2e7e | 33.57 | 312638 | 2.19 6.20 

30 | 32.52 | 18.14 | 2.41 3.96 | 34.69 | 11.86 | 1.94 8.04 

35 | 39.40 | 23.99 | 2.05 5.59 | 39.47 | 11.01 | 1.65] 10.12 

4.0 715 | 16.63 S29 1 3.77 1.36 | 24.86 | 11.85 | 2.99 3275 

20 | 20.39 7290 | 3.43 2.1C | 26.67 | 12.10 | 2.73 546 

25 | 24.84 | 12.77 | 3.09 3-13 | 28.85] 12.23 | 2.46 7.57 

30 | 30.05 | 17.83 | 2.73 4.51 | 31.66] 12.17 | 2.18] 10.01 

35 | 33.10 | 22.81 | 2.76 8.53 | 35.61 | 12.19 | 2.21] 19.11 

4O | 40.16 | 28.40 | 2.27] 11.96 | 36.36 |] 11.60 | 1.82] 23.32 

7.0 | 10 8.98 1.20 | 6.78 1.23 | 15.45 8.80 | 5.31 4.48 

13 § 32.22 5.48 | 6.02 2.40 | 17.10 9.52 | 4.73 8.38 

20 | 16.09 9.71 | 5.29 4.23 | 19.06 | 10.29 | 4.16] 14.00 

25 | 20.44 | 13.93 | 4.57 6.81 | 21.41 | 11.07 | 3.61] 21.0% 

30 | 25.26 | 18.26 | 3.91] 10.24 | 24.19 | 11.74% | 3.10] 28.98 

35 | 30.66 | 22.83 | 3.29] 14.70 | 27.53 | 12.17 | 2.62] 37.19 

4O | 36.98 | 27.84 | 2.70] 20.52 | 31.84% | 12.16 | 2.16] 45.35 

10.0 | 10 7-77 2.91 | 8.99 1.97 {| 12.01 7.09 | 7.01 770 

15 | 11.33 6.95 | 7.64 4.34 | 13.91 8.05 | 5.97] 16.32 

20 | 15.34 | 10.86 | 6.40 3.00 | 16.23 9.14 5.02] 28.56 

25 | 19.68 | 14.77 | 5.34] 13.06 | 18.90 | 10.23 | 4.20] %+3.36 

30 | 24.37 | 18.79 | 4.44] 19.70 | 21.92 | 11.21 | 3.51] 59.89 

35 | 29.57 | 23.07 | 3-68] 238.25 | 25.34 | 11.93 | 2.92] 76.84 

LO | 35.54 | 27.76 | 3.00] 39.25 | 29.51 | 12.24 | 2.39] 93.42 

45 | 42.81 | 33.20 | 2.37] 60.37 | 35.13 | 11.80 | 1.90 | 121.34 






































and lower limits on @ and in addition a value of .l/) be- 
low which no leading-edge separation is obtained. 


Lower Limit 


Across the reattachment shock the following rela- 
tions [Eqs. (157) and (155), NACA Rep. 1135] apply: 


_— MP&l(y — Ye + (y + D1 + 2G" — 1) 


M;? 
i? f ieee  — 1) 


(A.1) 
sin? 6 = [(y + l& + (vy — II (2yM1") (A.2) 


With these relations, Eq. (1), and Eq. (2) with? = 1, 
it is clearly possible, for specified 1/2, to solve for cor- 


responding values of .J/; and the turning angle 6, 
The lower limit on the compression corner angle, de- 
noted 8,,,;,, for leading-edge separation, may be visual- 
ized in the following manner. For a fixed 6 and de- 
creasing .l/,, the slope of the separation region dimin- 
ishes. In the limit, therefore, the angle 6, approaches 
zero and the separated region coincides with the com- 
pression-corner surface parallel to the free-stream; 
thus 


Bin = by (A.3) 


This implies that in the limit, instead of a shock a 
Hence 


Mach line emanates from the leading-edge 
tor Binin 
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Pressure ratio across reattachment shock 


MW M, A.4) 


Upper Limit 

An upper limit, denoted @,,,,,, is reached if the strength 
of the reattachment shock is maximum compatible 
for maximum flow-deflection 
(168) of reference 4 


with a weak shock 1.e 
across the second shock. Eq. 
gives the appropriate relation, where the prime de- 


notes values at maximum flow deflection 
l 
ty { M, j ,? 


Viiy + I(y + ICI’)! + Sy — ICM’)? + 16]; 
A.) 


sin? 6.’ xX }(y + 1M’)? — 4 4 


Since Eq. (A.5) must be satisfied simultaneously with 
Eq. (1), and Eq. (2) with ¢ = 1 and with Eqs. (A.1) 
and (A.2), a unique set of flow values across the second 
shock exists. Carrying out the computations by suc- 
cessive approximations on a desk calculator we obtain 


the limiting values: 


M,’ = -1.1800 &’ = 13184 
6.’ 3.424 M,' 0.9539 
0.’ 72.654 


External Supersonic Flow 


SEPARATION "1 


These flow conditions must prevail across the re- 
attachment shock for every case in which 38,,,, is 
achieved. With .1/,’ fixed, we can calculate the lead- 
ing-edge shock angle for various J) by use of Eq 


(132) of reference 4 


M,’‘) 
y + 1)?.V/,* sin? 6 — 4(.W)? sin? 6, — 1) x 
(y.Wo? sin? 6; + 1) 
[2y.Wo? sin? & — (y — 1)\[(y — 1)? sin? 6 + 2 


(A 6) 


By rearrangement of Eq. (A.6), a quadratic in sin? 6; is 
obtained. Having solved for 6, 6; may be obtained 
from Eq. (2). The upper limit of the compression 


corner angle is therefore 


versus .\/) is shown in Fig. 2. Its 
indicates that there 


A plot of By. 
intersection with the curve of 8,,,,,, 
is a minimum value of the free-stream Mach number 
(Mo) min = 1.27 (for which 8 5.15°) below which no 


separation is possible 
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Thus 6 is a direct measure of the convergence (+6) or 
divergence (—6) per cycle. The dynamic and aero- 
dynamic representations are in all respects adequate for 
nonsinusoidal motions, except for the structural damp- 
ing. In the presence of structural damping, this in- 
terpretation of 6 must be limited to nearly sinusoidal 
motion, 6 = 0 (g). This is consistent with the usual 
experimental technique for measuring the structural 
damping of a lightly damped system from the decay of 
the unforced oscillations. Further, 6 cannot be inter- 
preted as a fictitious uniform structural or viscous 
damping required to produce sinusoidal motion, as the 
usual flutter solutions can be interpreted. Obviously, 
for 6 = 0 the motion is sinusoidal and the critical speed 


is correctly determined for any value of g,,». 


0 the roots of the characteristic 
will be 


In general, for g 
equation |A* in Eq 
complex, and in the mathematical sense there will re- 


(1.26), x? in Eq. (2.9) 


sult 

w +(2(1 + 126) A.2) 
which would imply instability. However, on physical 
grounds, for positive structural damping to result in 
damped motion, only the positive or principal root is 
significant [see Eq. | (2.13)]. On this basis, the sign of 
the complex part of the frequency parameter (A* or 
x”) implies the sign of 6 in the absence of aerodynamic 
damping. This is the criterion implied in the discus- 
sion of section | and in the calculations of section 3 








Blunt-Body Heat Transfer at Hypersonic 
Speed and Low Reynolds Numbers’ 


ANTONIO FERRI,* VICTOR ZAKKAY,** ann LU TING*** 
Polytechnic Institute of Brooklyn 


Summary 


An analytical method for the determination of the effect of 


shock curvature on heat transfer in the region of the nose has 
been developed. It is shown that for practical body shape the 
viscous terms in the Navier-Stokes equations are not important 
in the region of the flow far from the wall, and the displacement 
thickness can be neglected. Then the flow can be approximately 
represented by an inviscid-flow solution having as boundary 
conditions the body shape, which is not affected by the Reynolds 
number, and by a boundary-layer type of flow near the wall, 
having appropriate boundary conditions. This approach per- 
mits us to determine the heat transfer in the region of the nose 
even at very low Reynolds numbers. 

Experimental results are presented. The experimental results 


agree with the values given by the analysis 


List of Symbols 


A = quantity defined in Eq. (37) or Eq. (55 

C = pp/peete = constant = 1 

é = percentage error in velocity or velocity derivatives 
E = percentage error in pressure 


| = function defined by Eq. (22) 
= h,/h,,, Eq. (24) 
value of function g at the wall = hy hy, 


ll 


“eos 
ll 


uy /Uy 

h = local stagnation enthalpy 

h,, = reference stagnation enthalpy, Eq. (39 

h,, = free-stream stagnation enthalpy 

m = mass flow 

AJ, = free-stream Mach number 

p = static pressure 

Pr = Prandtl number 

gq = heat-transfer rate 

du = heat-transfer rate at the wall 

r = radial distance from the axis of symmetry 

R= reference radius of the body; radius of sphere, or radius 
of maximum cross section of body 

j = radius of curvature of the body 

R, = Reynolds number referred to stagnation conditions 
[Eq. (1)] 

s = function defined by Eq. (19) 

« = velocity component tangential to the body surface 

u, = reference value of u for boundary-layer calculations 

uo = value of uw at the wall for the inviscid-flow [defined by 


Eq. (2)] 
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“4, = value of uw at the matching point y = y 
layer and inviscid flow 


of boundary 


i = velocity component along the y axis 


ll 


velocity coefficient defined by Eq. (2) 


I’. = free-stream velocity 

\ = curvilinear coordinate tangential to the body surface 

y = curvilinear coordinate normal to the body surface 

¥; = station where boundary conditions for the momentum 
equations of the boundary layer are applied 

\ = station where boundary conditions for the energy equa 
tions in the boundary laver are applied 

= ratio of specific heats 

I = function defined in Eq. (54) 

\ = detachment distance = Yshock — Vbody 

€ = (y¥ — 1)/y + 1)}1 4+ [2/(—y -— 1) Mo} 

n = coordinate of boundary layer defined by Eq. (18 

# = angular coordinate, ~ x/R near the stagnation point 

= = 7 

A © Goon: / Peron [Eq. (56 

a = viscosity 

ume = reference viscosity 

“mo = viscosity for stagnation conditions 

0 = density 

p = density for stagnation conditions 

o0 = variation of stagnation enthalpy defined by Eq. (45d 
or velocity coefficient defined by Eq. (2) 

T velocity coefficient defined by Eq. (2) 

y = stream function defined by Eq. (22) 

w® = vorticity coefficient defined by Eq. (2) 

6* = displacement thickness 


(1) Introduction 


— RECENT INTEREST in hypersonic lifting vehicles, 
re-entry space vehicles, hypersonic gliders, and 
hypersonic airplanes has generated increasing interest 
in the determination of heat transfer at Reynolds num- 
bers sufficiently low so that the concept of classical 
boundary layers does not directly apply, but still suffi- 
ciently high to permit the use of continuum theory at 
least in the region near the stagnation point. 

In this range of Reynolds numbers and for hypersonic 
geometrical configurations of practical interest having 
blunt noses, there exists large interference between 
the viscous flow near the wall and the vortical flow 
produced by the attached shock at the nose. The im- 
portance of such interaction and the necessity for con- 
sidering it has been discussed for the first time in refer- 
ence 1. Investigations for the zero pressure gradient 
case and small vorticity have been published in refer- 
ences 2-8. 

In references 9-12 a discussion of such problems 
and a method of analysis for the stagnation point have 
been presented. In references 13 and 14 a method for 
the analysis of this effect at the stagnation point, and 
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in the region downstream of the nose, has been outlined 
by the senior author. In this report the analytical 
work summarized in references 13 and 14 is presented 
in detail and analytically justified. Some experi- 
mental results are presented. 

The basic approach discussed here is to retain and 
modify the boundary layer concept for the low Reyn- 
The flow field behind the shock 


one near the wall where 


olds number range. 
is divided into three regions: 
the boundary-layer-type differential equations with 
appropriate boundary conditions still apply; another 
outside this layer, where the viscous and conduction 
effects due to velocity and temperature gradients have 
small effect on the heat transfer at the wall and can 
be neglected in the first approximation; and a third 
which is the shock region, where the viscous and con- 
duction effects can also be neglected. It is shown 
that this approximate representation of the flow field 
gives the heat transfer at the wall with sufficient ap- 
proximation, even at very low Reynolds numbers. 


(2) Discussion of the Problem 


Consider an axially symmetric or two-dimensional 
body having a blunt nose and moving at hypersonic 
speed. curvature of the shock, the 
flow in the nose region is rotational. The shock- 
induced vorticity is a function of the flight conditions 
and body shape. The vorticity produced by the 
viscous effects at the wall is a function of flight condi- 
tions and Reynolds numbers. For high Reynolds 
numbers and for a cold wall, the flow field outside the 
boundary layer can be obtained by using the body shape 
as boundary conditions and by assuming that the shock 
The boundary laver can 


Because of the 


is a physical discontinuity. 
be determined by neglecting the presence of vorticity 
and displacement thickness, and by using as external 
conditions the pressure and velocity distribution 
determined in the absence of the boundary layer. 
When the Reynolds number decreases, the representa- 
tion of the flow field as a shock front with very small 
thickness, followed by an inviscid flow and by a bound- 
ary layer flow, must be questioned. When the shock 


is considered as a physical discontinuity, the shock 
shape defines the conditions of the flow behind the shock 
by means of the shock conservation equations. These 
equations state that the velocity component tangent 
to the shock and the stagnation enthalpy remain con- 
stant throughout the shock. The boundary conditions 
required at the shock, for the flow behind the shock, 
are that the velocity at the boundary shall match the 
velocity given by the shock relations, and the stagnation 
enthalpy shall be equal to the stagnation enthalpy of 
the flow in front of the shock. No consideration is 
required for the derivatives of stagnation enthalpy and 
velocity components in normal direction to the shock 
When the Reynolds number decreases, and the 
shown (refer 


front. 
shock thickness increases, 
ence IS) that the existence of derivatives of velocity 


it has been 


and stagnation enthalpy in normal direction to the 
shock, produces appreciable variations in the tangential 
velocity component and stagnation enthalpy; and 
therefore, in order to define the flow properties behind 
the shock, velocity components, total enthalpy, plus 
the derivatives of velocity components, and _ total 
enthalpy in normal direction to the shock must all be 
specified. The effect of variation of velocity deriva- 
tives and total enthalpy behind the shock on the values 
of the velocity components and total enthalpy, can be 
related to the Reynolds number and to the values of 
derivatives. It has been shown that for practical 
applications this effect is of the order of 


1 Ry + 1) (y -— 1 
where 
R. = (poRYV hy.) / pe (1) 


is the Reynolds number related to stagnation conditions 
(reference 1S). This implies that the errors intro- 
duced in the flow field by using the shock conservation 
equations are smaller or on the order of 15/R,; and 
thus modifications of the classical shock equations, in 
order to determine heat transfer on the body, are re- 
quired only at very low Reynolds numbers when 15/R, 
—> |. 

Let us now consider the flow between the shock and 


the boundary laver. For boundary conditions of prac- 
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tical interest, it is of importance to analyze in detail 
the properties of the flow field in this region. The 
velocity components 1 and v, tangent and normal to 
the body surface, are functions of the body shape and 
of the free-stream conditions, and can be expressed as 


follows: 


u = uo(x)[1 + w(y/ R) + o(y?/R*)] +... (2) 
v = v(x) [(y/R) + r(y?/R2)) +... f 1“ 


where x and y are curvilinear coordinates along and 
normal to the body and & is a reference length for 
example, the radius of the body. For bodies of prac- 
tical shape, in the region of the nose the variation of 
velocity components is almost linear with y, and the 
quantities o and 7 are small. As an example, in Figs. 
la and 1b the distributions of u and v as functions of 
y/R for different values of x/R are presented for a 
spherical body of radius R at 7 = 20 and 200,000 ft 
altitude. The values of w, (Om/V..)/(Ox/R), o, and 
7 for the distributions of Fig. 1 are also given in the 
figure. The value of w is a function of the density 
ratio through the shock, and of the body shape—for 
example, for a spherical body at high Reynolds num- 
bers the quantity in the stagnation region mp is given 
approximately by the modified Newtonian distribution 
and can be expressed as 


to = V.V[(y — 1)/y]8 (3) 


> / “ ; 
where 0 = x/Rand I’. ~ V 2h,,. Then the quantity 
can be expressed approximately as 


we ={Viy/(y- I) - 110 - ©) 
[el — V8/3 e+ <a (4) 
where ¢€ = [(y — 1)/(y + 1)] X 
11 + [2/(y — 1))/M.?} 


because in the region of the nose the shock is close to a 
spherical shock. For this body the quantity v has 
the form 


{1 + [(A/R)r]lo = —V.[1 — (82/2)] (5) 


At low Reynolds numbers the viscosity and conductivity 
of the gas modify the flow field and therefore modify 
the boundary conditions outside the boundary layer 
and affect the value of the heat transfer at the wall. 
For the case of laminar boundary layer, the heat transfer 
varies with +/m and h,,, where 1 is the velocity and 
h,, the stagnation enthalpy external to the region con- 
sidered as boundary layer; the heat transfer is effected 
by the value of 0u/Oy, and slightly by the value 0,/Ow. 
At low Reynolds numbers the effect of 0u/Oy can be 
as highas 50 percent. In order to determine the effects 
of viscosity and conductivity in the heat transfer it is 
necessary to determine the effect of these quantities 
on the values of 1, hy,, (Ou/Oy);, and Ou,/Ox as func- 
tions of the Reynolds number. An error of the order 
of e in the quantity , due to incorrect consideration 
of the viscosity and conductivity effects, will introduce 
errors of the order of 1/2e in the value of the heat trans- 
fer. It will be shown later that an error of e in the 
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value of Ou;/Oy produces an error much smaller than « 
in the value of the heat transfer. An error of the order 
of ec in the value of /,, will introduce an error of the same 
order in the heat transfer. The order of magnitude of 
the errors introduced by neglecting conductivity and 
viscosity in the quantities 1, h,,, (Ou/Oy),, and (0u/Ox), 
in the region outside the boundary layer, can be ob- 
tained by determining the order of magnitude of the 
differences between the values given by the Navier- 
Stokes and Fourier type of equations and the values 
given by the momentum and energy equation in ab- 
sence of viscosity and conductivity. The leading 
terms containing viscosity in the Navier-Stokes equa- 
tions can be expressed as (reference 20): 


) 4a — r) + (w = 3) [(Or0 Ov) + 11g cot |} Mo R? 
(6a) 


(Ovo /O#) | o/R? (6b) 


I2u9a0 + 2uow — 3p — 
The first term effects the quantity 0p/0r; therefore the 
error / on the value of the pressure at the surface of the 
body introduced by neglecting this term is of the order 


of 


Ano 


R | se + T) + 


( 3) (oe i» at 
w® e oe Up COT t 


where A is the detachment distance and A/R ~ 
(y — 1)/(y + 1). The percentage error ¢ in velocity 
u is of the order of 1/2#, therefore the percentage 
error in velocity is of the order of 


; ] Ay 7 P 2’ Olp 
é= > RY, 4oo(1 + 7) + (w — 3) ay + mM cot & 


p:' —~p. = Ep ~ 


but 
Ps ™ PaVn*, po/p. ~ (y — 1)/(7 + 2), 


V. = V2h,, and R, = (pV hyR)/po 
therefore 


4v 


wal a+ 0 + —-x 
¢~ SeaJaivy."** 


I 


(Se . - ais 
v Ug COt 2 (4) 


This term is small—for example, for a spherical body 


the quantity e becomes 
e~ {V[(y—1)/y](w — 3) — 2(1 + 7)(w/V.)} 
R.vWf2 (8) 
If y is of the order of 1.15 to 1.25, then 
e~ 10/R, 
The term of Eq. (6b) effects the quantity (1/R) (Op 
Ov) which is of the order of pv(Ou/OR), p(udu/ROdI), 


(uv/R)p; therefore the percentage error e in u, given 
by neglecting this term, is of the order of 
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Fic. 2. Displacement thickness as a function of stagnation 
Reynolds number at a station 45° of a spherical body (cold 
wall hypersonic flow ) 


] Ov | 
e= 260 + 2w — 3 ++ 
2 Ov Io 
Mh 


pol(Ouo/O9) [1 + w(y1/R)} 


(9) 


where y; is the station corresponding to the outer edge 
of the boundary layer where the boundary conditions 
are imposed. The value of e given by this expression is 
also small. For example, for a spherical body at hyper- 
sonic speed, and low Reynolds numbers y; > 1/2 where 
A is the shock detachment distance, this distance A is 


of the order of 


A/R = e{1 — V8/3e + 4e] (10) 
therefore 
cre (20 + 20 - 3 + ) x 
2 \ Wy 1 
V y¥/[2(y — 1)] Lan 


Vil + w/2(y¥ -— D/(y + 1) 2! 
Then for values of y of the order of 1.15 to 1.25, 
e < 10/R, (12) 


In a similar way it can be proved that the error in the 
value of the stagnation enthalpy h, outside the bound- 
ary layer is also small. Express the stagnation en- 
thalpy in presence of viscosity and conductivity as 
h,’ = h,,(1 + e); then for finite values of the coordinate 
3 the order of magnitude of the error ¢e is given by 


Oc R- By” 
ov h,, pu R* 
, (1 — P ) (Ov Ov 
| Pr uy SJ 


where Pr is the Prandtl number. The value of e given 
by this expression is small—-for example, for a spherical 
body, where uo and v% are expressed as in Eqs. (3) and 


(5), the value of e is given by 


ve VAy — 1)/¥ 
e —= 


ss 9 R | ee + i) = 


Pr—1 [oy | 
———} — (13 
( rv ) we= 1 } 
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If Pr = 0.75, y = 1.22, and w is given by Eq. (4), then 

e~ 25 8 /R, (14) 
The error is of the order of e < 15/K, for 3 between the 
stagnation point and the sonic point. 

Near the stagnation point, where # — 0 and v > 0, 
the leading term in the energy equation is yu O*h,/Oy’, 
and the error due to conductivity is independent of the 
Reynolds number and is of the order of 


po (2(Ouo Ow’)? Vo" | oe 1.5A° 
poA~*h,, R? R? 


e¢~ = 3 (15) 
When v is different from zero and 3 — 0, the leading 
term in the energy equation becomes the term pv(0h,/dy), 
and the error becomes of the order of 


Ouo(OlMo/ Od)” 
e~ : (16) 
pol’ h,,.R ; 


2 


for spherical body e ~ 1/K,. 

From these considerations it appears that, for bodies 
of practical interest, the effect of viscosity and con- 
ductivity in the flow outside the boundary layer will 
modify the heat transfer at the wall between the stag- 
nation point and the sonic point of quantities of the 
order of 15/R, or less, and therefore the percentage 
errors introduced in the evaluation of heat transfer 
if the influence of viscosity and conductivity is neg- 
lected in the external flow is of the order of 5 percent or 
smaller for Reynolds numbers as low as 200 to 300. 
This value of the Reynolds number is of the same order 
of magnitude as the Reynolds number below which the 
use of continuum theory is open to question. From 
this analysis it appears that from the point of view of 
determination of heat transfer, the introduction of the 
Navier-Stokes-Fourier type of equations in the flow 
between the shock and the boundary layer does not 
seem useful, unless all the assumptions introduced in 
the calculations have consistently the same required 
precision. The limit of continuum theory can be ex- 
pressed by the condition® 


,/A <1 


where A, is the mean free path behind the shock, now 
in the region of the stagnation point 


\, M “or i 
~ - a P (17) 
A pe VOR (yy —-DV(y. -— DR 
If the condition \,/A ~ 1/10 is assumed as a limit of 
applicability of continuum theory, X, is of the order 

of 160 to 250 for values of y between 1.15 and 1.25. 
Now let us consider the boundary layer near the wall 
for conditions of decreasing Reynolds numbers. When 
the Reynolds number decreases, the shear stress at the 
wall decreases and at low Reynolds numbers becomes of 
the same order as the shear produced by the shock- 
induced vorticity (reference 13). Therefore, this effect 
must be included in the boundary conditions for the 
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boundary layer, and new boundary layer analysis must 
be developed. Because of the decrease in Reynolds 
number, displacement boundary-layer thickness and 
boundary-layer thickness change. The displacement 
thickness 6* for cold wall is very small, and the changes 
with Reynolds numbers are not important. As a con 
sequence, in this approximation the shock shape is not 
changed and the pressure distribution along the body 
surface is not affected by the change in Reynolds num- 
ber. In Fig. 2 the displacement thickness as a func- 
tion of the Reynolds number for a sphere at a station 
15° from the stagnation point is shown. Newtonian 
flow is assumed for the pressure distribution and the 
wall temperature is assumed to be zero. 

When the Reynolds number decreases, the boundary 
layer thickness increases. The increase of thickness of 
the layer where viscous effects are considered intro- 
duces two different limitations. In order to use 
boundary layer approximations it is necessary that (1) 
the ratio 6/r of the boundary layer thickness (6) to the 
radius of curvature of the body (7), be much smaller 
than 1; and (2) the boundary layer thickness be smaller 
than the detachment distance A. The ratio A/F, where 
R is the radius of the maximum cross section of the 
body, is small. For a spherical body, 
AR ~ (y — 1)/(y + 1) 


A/R 0.1 for y= 1.22, 
and A/R ~ 0.07 for y =1.15 


The condition 6 < Ais therefore representative for both 
limitations. In the analysis presented in the following 
section, the flow field behind the shock is divided into 
two regions: one near the wall where boundary layer 
approximation is used, and another near the shock where 
conductivity and viscosity are neglected. Velocity, 
velocity derivatives, enthalpy, and derivative of en- 
thalpy are continuous at the matching of the two re- 
gions. The station y, where the matching of the two 
velocity profiles occurs is a function of Reynolds 
number of the velocity distribution along the body, and 
of the value of y. In Fig. 3, the values of the Reynolds 
numbers corresponding to y, = A and y = 3/44, as 
functions of y and determined with the method pre- 
sented later, are shown. The body considered is a 
sphere. If the body is different from a sphere and the 
radius of curvature at the nose increases, then the value 
of A/R is larger than (y — 1)/(y + 1), and therefore 
the value of KX, corresponding to y, = A or y = */4A 
decreases. The curve shows that for values of y of 
the order of 1.25 to 1.15, which are realistic values for 
hypersonic conditions, the &, where y, < A is of the 
order of 100. This implies that the boundary-layer 
approximation can also be used with sufficient accuracy 
at very low Reynolds numbers close to the limits of the 
continuum theory. From these considerations it ap- 
pears that it is possible to analyze the flow behind a 
detached shock even at very low Reynolds numbers by 
introducing simplified assumptions for the shock and 
the vortical layer without great loss of precision. In 
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this approximate scheme of analysis the shock is repr 

sented as a region across which the stagnation enthalpy 
and tangential velocity component are preserved 
The flow in the region behind the shock has a constant 
stagnation enthalpy and zero gradient of stagnation 
In this region 
viscosity and conductivity are neglected. The deter- 
mination of this region of the flow field is independent 


enthalpy normal to the shock surface. 


of the Reynolds number and is performed by deter 
mining the flow field for a body shape equal to the actual 
body with the same methods used at high Reynolds 
numbers. Near the wall a boundary-layer type of 
analysis is used with boundary conditions appropriati 
to the low Reynolds number case. 

Some different types of analyses have been presented 
for the stagnation point on a cylinder or a sphere in 
references 9, 10, and 12. Here four regimes have been 
defined for Reynolds number ranges corresponding to 
conditions for which the fluid can be considered to be 
continuous. Several assumptions have been intro 
duced in the analysis for the shock shape and flow dis- 
tribution. Because of these assumptions on the bound 
ary conditions it is not clear whether the use of dif- 
ferential equations that retain higher order terms will 
really improve the precision of the results. 


(3) Boundary Layer in the Presence of Shear 
Flow 
A. Differential Equations 
Following the method used by Lees in reference 15, 
consider two coordinates 7 and s, defined by 


~~ 
ell, 
n = £ 1/2 | ry" 4 dy (1S) 
toa? JO Pp, 
*X 
s = | Ph U,ry* dx (19) 
/J0 


where all the quantities with the subscript e are inde- 
pendent of y and are functions only of x and will be 
defined later. The quantity & is equal to 1 for the 
axially symmetric case and is zero for the two-dimen- 
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Fic. 3. Reynolds number corresponding to the conditions 
yi = Aand y, = °/, A on a spherical body as a function of - 


(hypersonic flow). 
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sional case; the quantity 7) is the radial distance from 
the axis of symmetry. 

If the subscript 1 indicates the conditions immedi- 
itely outside the region considered as the boundary 
layer, then 


pit; du,;/dx = —dp,/dx (20) 
\ stream function is introduced, defined by 
purs’ = Oy/dy (21a) 
puro’ = —Oy/dx (21b) 
and the function y is expressed as 
(Sn) = f(n)(2s)¥/? (22) 


The momentum equation for a boundary layer flow 
then becomes 


(cf’’)’ +H" + (2s/u,)0u,/Os[pi/p( fr’)? — 


(f’)?] = 0 (23) 

where ¢c = pu/p.u,. The velocity u becomes 

u = u,f'(n) 
and the quantity /,’ corresponds to #/1,. 

Let us define 

g(n) = h,/h., (24) 
The energy equation for Prandtl number | and c = | 
then becomes 

g" + fe = 0 (25) 


The momentum Eq. (2%) differs from the equivalent 
equation of reference 15 in the last term of the equation. 
This term has the same order of magnitude in both 
cases and similar variation with respect to 7 (reference 
20). It has been shown in reference 15 that the con- 
tribution of the last term is very small and can be 
neglected without introducing large errors in the re- 
sults. Because the terms are of the same order, the 
last term of Eq. (23) can be neglected here on the same 
ground. 

Eq. (23) then becomes the well-known equation 


fr’ +f'f =0 (26) 


for which solutions are available. 


B. Boundary Conditions for the Momentum Equations 
Consider now a boundary profile in the presence of 

shear flow. At a given value of x or s, the boundary 
conditions at the wall are the standard boundary condi- 
tions of zero velocity, and at » = 0, they must be 

f(0) = 0, f’(0) =0 
The velocity distribution near the wall, given by the 
inviscid flow, can be expressed as 

u = Uo[l + w(y/R)] (27) 


where uw is the velocity at the wall and is a function 
of x, and wu is the velocity at a distance y from the wall 
in the inviscid flow. 
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At a given station x the matching of the boundary 
layer profile with the inviscid flow will occur at a sta- 
corresponding to a distance 
At this sta- 
tion the velocity u and velocity gradient Ou /Oy, given 
by the boundary profile and by the inviscid flow, must 
be the same, and the the 
boundary layer profile and between the stations 7 = 0 


tion 1 where uv = 1%, 
from the wall for the inviscid flow field. 


mass flow contained in 


and 7 = m in the inviscid flow, must be the same. Let 
uy) = Uo[1 + w(y/R)] = woG (28) 
because 
u/ue=fi', Uu, uo (fy! (29) 


In the boundary layer profile 


»yyrr 


( pil,*ro*) [(2s)} =I 


Ou/Oy = (Ou/On)(On/Oy) 


Therefore if the value of Ou /Oy must be the same in the 
boundary layer and in the inviscid flow, 


, 


( pi2t,?79") [(2s)! *\hi i= 


> 


uUqgw/ R (30) 


In the inviscid flow the variation of density from the 
value p; at the station y, to the value of pp at the station 
y = Ois negligible (u is small); therefore, the value of 
po can be used in place of p;, and from Eqs. (29) and 
(30), there results 

f/f’? = [w(2s)"/?]/ (ro Ruopo) (1/G?) (31) 


By definition, at the station y, it must be true that 


| e . lo + My 
y = (2s)!/2f, = ro’ pudy = ( : ) poViro’ 
J Z 


uo/2(G? — 1)(ro*¥Rpo)/w (32) 

Dividing Eq. (31) by Eq. (32), it results that 
1/G? = 1 — [(2ff'")/(fi’)?] (33) 

The combination of Eqs. (31) and (33) gives 
(ro*uopol)/[w(2s)"/27] = (fi'2/fi’”) — 2h (34) 


Eq. (34) defines the value of the station 7, where the 
boundary profile defined by Eq. (26) matches the ex- 
ternal flow defined by Eq. (27). 

At each station x the quantity in the left side of Eq. 
(34) can be evaluated; when the inviscid flow is known 
the quantity in the right side is a known function of », 
obtained from tabulated values of the solution of Eq. 
(26); therefore, the value of » m can be obtained; 
then Eq. (33) defines the value of #, and Eq. (29) the 
value of u,. 

In the region of the stagnation point, the quantity x 
can be expressed as 


x = KO (35) 
where # is the angular distance from the stagnation 
point. Then 

ro = Rsin d => Kd (36) 
The quantity can be expressed as 
Uy = AVh,, 3 (37) 
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where 


Poy = Pebe 


» 


Then by introducing Eq. (33), expression (34) becomes 
near the stagnation point, 


IVR.V(k + 1)A]/o 
(fi’ —_ 2h f°") . I ( f,’2 hy 


? 


) — Zi} (38) 
where X, is the Reynolds number defined in Eq. (1) 


C. Boundary Conditions for the Energy Equations 

The energy equation in the boundary layer for the 
case of very cold wall can be expressed in the form of 
Eq. (25) which permits a solution 


h./k., = ¢ = f' 29) 
The boundary conditions at the wall require 
£(0) = Su (40, 


At the matching point 7 = m» of the viscous and con 
ducting layer with the inviscid flow, two boundary 
conditions are required: 


2(m) = h,,/hs,, and g'(m) = 0 (41) 


At the station »; where the inviscid and viscous velocity 


wrofiles are matched, g(n,) f'(m,); therefore 
I £\n m1), 


(Og/Oy),; = 1/u,(Ou/Ov), # O 


If the boundary conditions for the energy equation 
are applied at the station » = m, a discontinuity is 
introduced in the value of 0h,/0y which is not physical 
This discontinuity is on the same order as the vorticity 
term in the velocity distribution, and thus cannot 
be accepted. The boundary conditions for the energy 
equation cannot be imposed at a station » = m, but a 
region above n; must be considered as part of the bound- 
ary layer where conductivity is important. A general 
solution of Eq. (25) has the form 
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an sun 
y = c| dn Jo + gu 


, (42) 


So’ 
, ’ 0 


Y = Ce 
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This solution indicates that the quantity g’ decreases 
rapidly with increasing values of » when » < m. This 
implies that the conductivity effects must be con 
sidered only in a small layer outside m between a sta- 
tion m and yn defined by the conditions of Eqs. (41) 
expressed in the form 

2(m2) = hy, /h., and g’(ns) <K (43) 
where ¢ is a small value. 

The presence of heat conduction in a small layer \ 
between station y,; where the velocity profile is matched, 
and a station y. where the energy equation is matched, 
changes the density distribution, and therefore affects 
The effect on the 


pressure can be analyzed in terms of changes of dis 


velocity and pressure distribution. 
placement thickness. This effect is small and can be 
neglected (reference 20). If the pressure distribution 
remains unchanged in the approximation accepted, 
ind uw’, uw’, and p’ are the quantities for the flow with 
heat conduction, then 


Ou 9) ( | , , on’ a) ( oe) (44) 
l = ~pit = 
mm Ox Ov * Ov ee” Ox ov ‘ Oy 


Assume 


u’ = u(1 + e) oI (w/R)yv\(1 + e) (45a) 
Ou’ Ov (1tywo/R)(1 + 7) (45b) 

O71’ /Oy? T (45c) 

h, = h,,\1 — a) (45d) 


Then 


(0a/OY) y=y, — (yw) /(u,R) (46) 


and at vy = \ it must be 


(0a/Oy),_,. ~ O 


Express 
o (wity/2Ru,)(v2 — v)?/r (47) 
where A = y. — \. 
Then in the approximation accepted, it is 
a|pu(Ou/Oxr) + (tw/ R)w/A) SS — 2epu(Ou/Ox) + ry 
(48) 


Near the stagnation point at low Reynolds number \/R 
can be of the order of 0.5; w can be of the order of 20. 
Then 


o wig 


2 tn, R 


lly ? rN R, ¢ — 
r = ( ) Se gly — 1) (4g) 
u.J R4w \ Y 


Then the effect on heat transfer introduced by neglect- 
ing viscosity in the layer of thickness \ is small and, for 





Vv 


V 


Vv 





464 


BLUNT-BODY HEAT TRANSFER 


tegration for the region between my, and y can be per 
Then at the station m, 
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ct- wall and g; and g;’ at 7; can be determined. Then the Fic. 6. Distribution of gw vort./qu vor 0 along the spherical 
£1 o_o : ermined ty . ' body at h = 300,000 ft and Wa = 20and R, 1,470 for R spher. 
A numerical in- 2 ft 


for constants in Eqs. (42) are known. 
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Fic. 7. Photograph of models. 


or with better approximation on the basis of the analysis 
of reference 16. The ratio A Qwvorticity /Qwvorticity—o 18 


given by 


/ ~~“ 

A = qQvvort./Quvort.-0 = V te/Ug(he,/ he) (56) 

In Fig. 4, three values of h,,/h,, have been considered, 
hy,/hs, = 0, 0.84, and 0.46. The value of /,, has very 


small effect on the quantity A; however, it has a large 
effect on the value of , where the matching between 
external flow and viscous flow occurs. The values of 
m as functions of IT are also given in Fig. 4 for the three 
values of /,/h,, considered. 

The quantities w and .1 are functions of the value of 
y and of the body shape. The relation between A and 
I given by Fig. 4 permits the variation of A as a func- 
tion of the Reynolds number for a given body shape 
and value of y to be obtained. 

In Fig. 5 results of such calculation are shown for the 
condition « = (y — 1) (vy + 1) = 0.1 and for a spheri- 
cal body. In this case, the quantities w and .1 are 
given by Eqs. (3) and (4). 
Reynolds numbers y./k < AR, therefore, the condi- 
tions expressed by Eqs. (43) are used on boundary con- 
In the lower range 


For large values of the 


ditions for the energy equations. 
of Reynolds numbers the value of y./ obtained from 
these boundary conditions is larger than A, 2; therefore 
the boundary conditions given by Eq. (53) have been 
imposed at the station Ak. The lowest Reynolds 
number of the curve corresponding to the condition 
yi/R = A/K is also indicated in the figure. For com- 
parison the results obtained for the same conditions by 
methods of references 9-12, and presented in references 
10-12, are also given. The trend given by these types 
of analyses is the same; however, quantitatively the re- 
sults are different. The results for the viscous layer 
assumption require the solutions of a system of Navier- 
Stokes-Fourier type of equations. 

The present method can be used in the regions far 
away from the stagnation point, with the same limita- 
tions of the method of reference 15. In Fig. 6, the dis- 
tribution of A as a function of the angular coordinate 3 
is given for a spherical body at the altitude of 300,000 
ft and M. = 20. The value of y behind the shock is 
1.14 and the Reynolds number Xk, = 1,470. As indi- 
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cated by the figure, the effect of vorticity on heat trans 
fer decreases gradually moving away from the stagna 
tion point. 


(4) Experimental Investigation 


In order to determine experimentally the effect of en 
tropy gradient on the heat transfer in the stagnation 
region, an experimental program has been initiated 
in the hypersonic wind tunnel of the Polytechnic Insti 
tute of Brooklyn. 

The experiments presented here have been per 
formed in the Mach 8 wind tunnel for two different 
stagnation enthalpy conditions. Because of the rela- 
tively small effects due to entropy gradients expected 
for the conditions of tests, a special technique has been 
developed in order to perform the experiments and to 
analyze the data. 


A. Experimental Technique 


Standard heat-transfer measurements do not permit 
very accurate results to be obtained because of com 
bined errors in the measurements of the stagnation 
conditions and in the heat transfer to the model due 
to model instrumentation. In the present investiga- 
tion heat-transfer data for a large range of Reynolds 
numbers are required. A variation of Reynolds num- 
ber of the order of 20 to 30 can be obtained in a given 
model by changing the stagnation conditions of the 
wind tunnel; therefore, for this range only, errors in 
free-stream flow properties are introduced. However, 
in order to obtain variations of Reynolds number of the 
order of 10%, several models must be used in the tests 
which introduce additional instrumentation errors. 

The change of tunnel Reynolds number affects 
slightly the free-stream Mach number and makes the 
analysis of the data more difficult. In order to elimi- 
nate or reduce these sources of error, the experiments 
have been conducted in the following way. In each 
test several models of different sizes have been placed 
on the test section. In this way a range of Reynolds 
numbers varying in the ratio from 1 to 20 has been 
covered in each given test. The heat transfer at 
equivalent stations in each model has been determined 
simultaneously. The models initially had exactly the 
same temperature, and the heat transfer has been ob- 
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c 


trans- 
tagna 


of en 
lation 
tiated 
Insti- 


per- 
‘erent 
rela- 
ected 
been 
nd to 


-rmut 
com 
ation 
due 
tiga- 
10lds 
1um- 
iven 
the 
‘S$ in 
ver, 
' the 
ests 


ects 
the 
imi- 
ents 
‘ach 
iced 
Ids 
een 
at 
ned 
the 
Ob- 





BLUNT-BODY 


tained from the measurement of variation of tempera- 
ture with time at the starting of the tests where the 
temperature is still uniform. The starting time of the 
tunnel is very small because the tunnel initially is 
evacuated, and a quick-opening valve is used which 
plugs the throat of the tunnel. 

Initially, tests at very high Reynolds numbers, where 
no entropy effects are expected, have been performed 
on models of large size. The value of the measured 
heat transfer, multiplied by ~/R where RF is the radius 
of the body, is constant in the range of high Reynolds 
number, and does not change from one model size to 
the other, confirming the well-known result that in 
the absence of vorticity, the heat transfer is inversely 
proportional to the square root of the radius of the 


body. Then several models having different sizes have 
been tested simultaneously at different Reynolds num- 
bers. The results obtained in each test have been 


inalyzed by plotting the ratio of the quantity gVR, 
where g is the measured heat-transfer rate for a given 
model, to the value of g\/R measured for the largest 
moel used in the same test. The ratio is unaffected 
by small errors in the free-stream Mach number and 
The effect of 
remains constant in 


due to 
instrumentation different 
This possible source of error has been eliminated 


stagnation conditions errors 
model 
tests 
by performing tests at high Reynolds numbers. In 
the larger models these tests permit Reynolds numbers 
to be reached where the value of gVR is constant; 
model instruments can be 
Then, by perform- 


therefore, a calibration of 
performed for the larger models. 
ing tests at lower Reynolds numbers with such models, 
the effects of entropy gradients have been obtained 
experimentally in the medium range of Reynolds num- 
bers. From these experimental results, possible in- 
strumentation errors on models having smaller radii 
have been determined. The tests have been repeated 
by interchanging the location of the different models 
used for the tests in order to eliminate possible errors 
due to variations of tunnel calibration. The models 
used in the present investigation have spherical noses. 
A photograph of some of the models used in the tests 


is shown in Fig. 7. 


B. Experimental Results 


The results of the tests are shown in Figs. 8, 9, and 
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on a spherical body at 7, = 2,300°R, M. = 8 
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10. The data in Figs. S and 9 give the effect of en 


The 
data of Fig. 8 correspond to a stagnation temperature 
of 1,600°R. 


tropy on the heat transier at the stagnation point 


The figure gives the ratio qwvort. Gwy 


as a function of Reynolds number. The quantity 
Gwvorr. iS the measured quantity, while the quantity 
Jpvort has been determined by reducing the quan 


tity measured in the largest model by the square root 
of the ratio of the radii and, if required, by the ratio 
Gwvort. Qwvort 
of the Reynolds nunber of the largest model in the 


given test. 


» determined experimentally at the value 


The second correction is very small and 
is necessary only in the lowest range of Reynolds num 
bers, where the Reynolds number of the largest model 
is below the value where vorticity effects can be 
neglected. 

In Fig. 9 similar results are shown for a stagnation 
temperature of 2,300°R. The 
temperature affects the values of the vorticity and 
shock detachment and, therefore, affects the value of 


variation of y with 


the ratio gwvort. Gevo 

In the figures the values of X, corresponding to the 
condition yo A and y, A are also given. The ex 
perimental results presented are all in the range of 


rt () 


Reynolds numbers where yo < A. 

The results have been compared with the results of 
the present analysis and with the analysis for vorticity 
interaction presented in references 9 and 12. The 
agreement of the experimental results with the present 
analysis is very good, while the analyses of references 
9 and 12 underestimates the effects of vorticity inter 
action. 

In Fig. 10 the ratio gwvort./Qwvort 
tions around a spherical body at two values of Reyn 
For 


o for different sta 


olds number is_ presented. comparison, the 
values obtained with the present analysis are also pre- 
sented. Again the agreement between experimental 


and theoretical data is satisfactory. 


(5) Conclusions 


An analytical method for the determination of the 
effect of shock curvature and entropy gradient on heat 
transfer on the nose region of blunt bodies at hyper- 
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The Effect of Sweep Angle on Hypersonic 
Flow Over Blunt Wings' 


SEYMOUR M. BOGDONOFF* anp IRWIN E. VAS** 


Princeton Unwersity 


Summary 


A series of tests were carried out in the Princeton University 
helium hypersonic wind tunnel on blunt two-dimensional wings 
at zero angle of attack with sweep angles up to 70° at Mach 
numbers from 7 to 18. The leading-edge Reynolds number varied 
from 3,000 to 25,000 
compared with the simple summation of the theoretical inviscid 
effect blast 
number) added to the viscous effect (calculated as if no sweep 


The measured pressure distributions were 


(based on wave theory using che normal Mach 


were present). For the unswept wing, the slope of the pressure 
decay was reasonably well predicted by the theoretical calcula- 
tions. 
the pressure distribution due to changes in leading-edge Reyn- 


The viscous theory reasonably predicted the variation in 
olds number. By subtracting the theoretical viscous effects, an 
inviscid Mach number dependence of the 2.2 power was found 
as compared to the value of 2.0 predicted by the inviscid theory. 
The same approach for the swept wing did not give consistently 
satisfactory results. Deviations above and below the calculated 
value by as much as 40-50 percent were measured and there 
seemed to be no systematic variation with either Mach number or 
Reynolds number. At a constant high Reynolds number, it was 
found that the pressure distribution varied with the distance 
along the wing with an exponent between about —0.53 and —0.58 
except for a rather sharp decrease which occurred for the 70° 
sweep case. The pressure at a given station for a fixed Mach 
number and given leading edge thickness varied as the cosine of 


the sweep angle to the 1.1 power as compared to the 1.3 power 
predicted from general geometrical considerations 
Symbols 
Cc = proportionality constant used in relationship w/u = 
Cc T, 7; 
Cy = leading-edge drag coefficient 
MJ, = free-stream Mach number 
p = measured pressure on the plate 
Pi = free-stream static pressure at the leading edge of the 
model 
p,' = free-stream static pressure at an orifice if no model 
were present 
Ap =p- pi’ 
Re, = free-stream Reynolds number per inch 
Rey, = Reynolds number based on the leading-edge thick- 


ness and free-stream conditions 
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Rey,, = Reynolds number based on the distance, 
stream conditions 
= distance from the center line of the leading edge along 
the model surface 


t = leading-edge thickness 
1\ = free-stream temperature at the leading edge 
7 = adiabatic wall temperature 
= ratio of the specific heats 
uy = sweep angle 
x = hypersonic viscous parameter 


Introduction 


Sian YEARS ago, the authors started a funda- 
mental study of hypersonic wings in an attempt to 
evaluate methods of extrapolating from low Mach 
number results to suborbital and orbital Mach num- 
bers.'~° The original work was particularly concerned 
with unswept wings; but some exploratory studies on 
delta wings were carried out in 1956 and 1957. Con- 
siderable difficulty arose in trying to analyze these 
results, the important 
variables such as leading-edge and apex effects. 
problems dictated a more fundamental approach and 
more detailed studies have been undertaken 


particularly in separating 


These 


several 
recently.°’ The present paper is a summary of these 
two most recent studies and is particularly concerned 
with the blunt two-dimensional wing, both unswept 
and swept, at zero angle of attack. 

The basic problem under investigation has been 
understood, in a general way, for some time. For 
hypersonic flow over an unswept plate, it has been 
found that there is a very strong leading-edge effect, 
essentially inviscid, which has been treated theoretically 
in the last few years by a “‘blast wave” analogy. This 
leading-edge effect is found to be quite strong at high 
Mach numbers and extends for many leading-edge 
dimensions downstream. The effect appears to persist 
down to leading-edge dimensions which are of the order 
of a few mean free paths. The effect is purely a geo- 
metrical one and depends only on the ratio of the 
linear dimension to the leading edge thickness. On the 
basis of the blast wave analogy, the effect is a function 
of M?.®° 
the boundary-layer growth, which essentially modifies 
Such theories have been available 

They depend on the Reynolds 


The viscous hypersonic problem depends on 


the body boundary. 

for almost ten years. 
number based on the linear dimension and vary as a 
function of .4*." Unfortunately, because of limited 
experimental data, it is difficult to evaluate these 
theories. The only proposal to handle the real problem, 
where there are both of the above-mentioned effects, is 
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Fic. 1. The Princeton University Helium Hypersonic Wind 


Tunnel 


to simply add the viscous and the inviscid predictions 
from the two theoretical treatments. '!)!* 

For the case of swept wings, the situation is even less 
satisfactory. The only proposed method of calculation 
is an extension of the above ideas of the unswept plate; 
that is, simply adding the viscous and inviscid effects. 
However, the effect of sweep angle must be included. It 
has been proposed that the blast wave analogy (inviscid 
effect) should be applied on the basis of the normal 
Mach number, the usual type of calculation at super- 
sonic speeds. The effect of leading-edge sweep on the 
viscous calculations is unknown and, to date, has been 
simply ignored. Such calculations have been made as 
if the plate were unswept; that is, the free stream Mach 
number and the streamwise dimensions have been used. 
Although this technique is easy to apply, it obviously 
has little to recommend it from the rigorous theoretical 
point of view. 

The present work is an attempt to examine, in detail 
and over a wide Mach number and sweep angle range, 
the two-dimensional blunt wing problem at zero angle of 
attack. The primary emphasis was on pressure distri- 
butions as a prelude to future heat transfer studies. 
The hope was to separate the viscous and the inviscid 
effects and to evaluate the calculation 
methods. The current program out at 
Mach numbers from 7 to 18, sweep angles from zero to 
70°, and for a range of leading-edge thicknesses to 
evaluate the leading-edge Reynolds number effects. 


proposed 
was carried 
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A wide range of exploratory studies of blunt wings, 
both swept and unswept, was presented in Ref. 1, which 
included results at angle of attack. The present report 
presents several major additions and modifications of 
these results. Previous data in the Mach number 17 to 
19 range were incorrect due to a combination of out- 
gasing and flow separation in the downstream section of 
the nozzle which was not obvious when the tests were 
carried out. The other preliminary sweep results of 
Ref. 1 at lower Mach numbers are correct for all com- 
pression angles and for low expansion angles; but the 
results at high expansion angles again suffer from 
problems of outgasing and are, therefore, not to be 
trusted. Further tests under these conditions will be 
necessary to completely validate the results for large 
expansion angles, particularly in the high Mach number 
regime. The present results have been extended down 
to a Mach number of 7 and the problems of outgasing 


and flow breakdown have been eliminated. 


Experimental Facilities and Models 


The test program was conducted in the Princeton 
University helium hypersonic wind tunnel'*® (Fig. 1 
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Fic. 2b. Axial Mach number distribution for the intermediate 


Mach number range, .V 11 to 14.5 


The tunnel is equipped with conical nozzles with re- 
placeable throat sections to produce Mach numbers 
from 7 to 9, 11 to 14.5, and 15.5to 19. Variation of the 
Mach number with a single throat was attained by 
varying the axial location of the model. The Mach 
number gradient along the axis of the conical nozzles 
was about 0.4 per in. for the low Mach number range, 
0.5 per in. for the intermediate Mach number range, and 
0.7 per in. at the high Mach number range. The 
centerline Mach number surveys for the three nozzles 
are shown in Fig. 2. The tests were carried out at 
stagnation pressures which varied from 150 psia at low 
Mach numbers to 1,600 psia at the highest Mach 
number. The free-stream Reynolds number varied 
from 0.5 X 10° per in. to 1.0 X 10® per in. Helium at 
room temperature was used as the test fluid and, since 
the run time was of the order of several minutes, the 
results presented are for the case of thermal equilibrium. 

The models were constructed of |, in. flat ground 
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steel, 2 in. wide, and about 35!) » 1n. long measured along 
the main orifice line. The thin leading and trailing 
edges were formed by grinding 10° wedges off the lower 
surface. The models were sting-mounted from the 
lower surface. The leading-edge thickness was varied 
by cutting off material normal to the test surface. The 
leading-edge thickness was measured, either under a 
microscope or optical comparator, and then formed to a 
hemicylindrical leading edge. The present study was 
made with leading-edge dimensions (diameters) of 
about 0.01, 0.02, and 0.03 in. The test program was 
carried out on plates with sweep angles of 0°, 30°, 45‘ 
52.5°, 60°, and 70° (Fig. 3). The main row of orifices 





Fic. 3. Photograph of the swept flat plates 
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(0.024 in. in diameter) was located in a region believed 
to be free from corner effects. Several static pressure 
orifices were located off the main row of orifices to check 
the two-dimensionality of the flow. The pressures 
were measured with an absolute silicone oil manometer 
using 20 microns as a base reference. The model was 
connected to the manometer with copper tubing to 
eliminate outgasing effects at the low measured pres 
sures. 

Although schlieren photographs were taken to check 
the general flow field, the usual ‘‘normal to the flow’ 
schlieren photograph gives no pertinent data for the 
swept leading edges. Previous studies'! included 
schlieren photographs taken parallel to the leading 
edge, but this work was not repeated. Several ex 
ploratory oil-trace studies were made to check the 
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general two-dimensional character of the flow and to 
obtain the boundary-layer flow direction. These 
studies were carried out by injecting a very low viscosity 
oil through one of the static pressure orifices located 
near the nose while the test was in progress. During 
this process, the oil streak was photographed. To 
facilitate the photography, the models were rotated 90° 
from their usual position so that the test surface was 
parallel to the windows. The models were sting- 
mounted so that the angle of incidence was always zero 
degrees. 

The complete study of the six sweep angles and three 
leading-edge thicknesses at Mach numbers of 7, 8, 11.6, 
12.7, 14, 17, and 18 are presented in detail in Refs. 6 
and 7 and only some typical results are presented herein. 


Experimental Results and Calculations 


The Unswept Wing 
A typical set of results for an unswept wing is shown 
in Fig. 4 in the form of log-log plots of A p/p,’ versus 
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s/t. From these data, it can be seen that the s/t 
exponent for the pressure decay along the plate is of the 
order of 0.55 to 0.62 with slightly higher values at a 
Mach number of 7. Other results for leading-edge 
thickness of 0.01 and 0.03 in. fall in the same range, 
considerably below the theoretical slope of 0.66 pre- 
dicted by the inviscid theory. 
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Fic. 6. Comparison between the viscous term subtracted from 
the total measured pressure, and the blast wave prediction 


The calculated pressure distributions were obtained 
as follows: The inviscid results, based on “‘blast wave”’ 
analogy, are only a function of the location of the orifice 
(rendered nondimensional by the leading-edge thick- 
ness) and the Mach number, and remain invariant 


versus s fas the leading-edge thickness is varied. 
Ap pi . C1M7(s t) 


where C, = 0.167 C,* 
and C,, is equal to 1.17 for a hemicyclindrical leading 
edge. For each value of leading edge thickness, the 
viscous effects were calculated by the theory of Lees- 
Probstein.'” 

The main parameter in their analysis is x defined as 


x = MiiVC V Re, 


When this parameter is small (weak interaction case) 


the pressures may be predicted by 
Ap/pr = C2 MPV C/V Re, , 
where 
“hy y [(0.865/M,") Ty 1, + 0.166 (y — 1)] 


There is, at the moment, no theoretical treatment to 
cover the real flow case where both viscous and inviscid 
effects occur. The only approach has been to add the 
two effects without considering any interaction 
effects. This approach has some validity if one of the 
terms is small compared to the other, but cannot be 
considered a valid theoretical approach; a “‘simple 
approximation’ is more descriptive. 

A comparison of these simple calculations with the 
experimental data is shown for Mach numbers 7, 
12.7, and 17 in Fig. 5. Since the pressure distribution is 
plotted versus s/f, there is a single inviscid curve for 
each Mach number, and a viscous curve for each 
different value of leading-edge Reynolds number. 
Within experimental accuracy, the experimental and 
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calculated pressure distributions check within the order 
of 10-15 percent, the experiments usually falling above 
the theoretical predictions. 

The viscous theory can be tested by checking the 
correlation of the different leading-edge Reynolds 
number results. The inviscid pressure distribution can 
be inferred by subtracting the theoretical viscous con- 
tribution from the measured values. This technique is 
shown in Fig. 6 for Mach numbers of 7, 12.7, and 17. If 
the viscous theory were exactly correct and the experi- 
ments perfect, the data for the three leading-edge 
thicknesses should fall on a single straight line for each 
Mach number. As seen in Fig. 6, the correlation is 
reasonable (scatter of +5 percent). Subtracting the 
theoretical viscous effects from the measured pressures 
gives good agreement with the blast wave prediction 
(shown dotted) for the lower Mach numbers. However, 
at the higher Mach numbers, there seems to be a dis- 
crepancy of about 10 percent, indicating that the J/? 
variation proposed by the blast wave theory is incorrect 
(assuming the calculated viscous effect 7s correct). 
Using the experimental results from Mach 7 to 18, the 
exponent of the Mach number was found to be 2.2 for 
the best match. 

A check of the consistency of the results over a wide 
Mach number range can be made by including, with the 
present results, the studies of Henderson and Johnston 
at Mach numbers of 17 to 24.!° For a constant leading- 
edge Reynolds number of about 7,000, the variation of 
the pressure at several fixed s/t stations is shown in Fig. 
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our previous results at Mach numbers from 11 to 18 
and Henderson’s from 17 to 24 gave slopes higher than 
that shown, the addition of the results at the lower 
Mach numbers and the comparison of the full set of 
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data available shows a linear Mach number de- 
pendence of 2.4. It should be noted that this is a 
combined viscous-inviscid dependence and the same 
value of 2.4 holds for Re;,, of 13,000. Subtracting the 
theoretical viscous effects, as above, results in an inviscid 
dependence of 2.2 over the full Mach number range of 7 


to 24.6 


The Swept Wing 

A typical set of pressure distributions obtained for a 
swept wing, A = 45°, over a range of Mach numbers, is 
presented in Fig. S for one leading-edge thickness. In 
general, the effect of Mach number on the slope of the 
pressure distributions is not large (—0.56 at M = 
7 to —0.52 at WM = 18). The level of the pressures, 
however, varies by an order of magnitude. 

With no exact theoretical solution for the unswept 
wing, the only proposed method of solution for swept 
wings has been an extension of the “simple approxima- 
tion” or addition solution.''!” It has been proposed 
that the inviscid pressure distribution is a function of 
the normal Mach number and distance normal to the 
leading edge. In the free-stream coordinate system, 
using the sweep angle, A, the inviscid pressure distribu- 
tion is then 

Ap pi = C\M,? (S t) a18 (COs A)4 é 
and C) is the same as in the unswept case. It has also 
been proposed that the viscous term be calculated as 
though no sweep were present. This term then only 
depends on the free-stream conditions and distance 
along the surface in the free-stream direction. 


Ap/pi = C2 MV C/V Re, 


Following the method used for zero sweep, a “‘simple 
approximation”’ of the pressure on the plate could be 
made by adding the two effects. 

A comparison of the calculated results with some of 
the data at various combinations of sweep angle, Mach 
number, and leading-edge thickness is shown in Figs. 9 
to 12. For a sweep angle of 30°, the data exhibit little 
of the viscous effect calculated and the actual data fall 
considerably above the predicted results for the high 
Reynolds number case where the calculated viscous 
effect is about 50 percent of the inviscid calculation. 
This curve is typical of results obtained at other Mach 
numbers for this sweep angle. For a sweep angle of 45° 
(Fig. 10), the low Mach number results show excellent 
agreement with the calculations, but as the Mach 
number is increased, this agreement disappears, partic- 
ularly for the low Reynolds number case where the 
calculations underestimate the actual pressures by a 
considerable amount. At 60° sweep for WZ = 7, Fig. 
11, the low Reynolds number calculations are in good 
agreement with the experiments (calculated viscous 
effects about equal to the inviscid effect). The high 
Reynolds number calculations fall considerably above 
the experimental results. The variation in pressure 
measured due to the viscous effect is considerably larger 
than predicted. For the same 60° sweep angle at a 
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Fic. 13a. Pressure distribution on a flat plate for varying A, 
1 = 7, Rei,, ~ 9,000 


Mach number of 12.7, the agreement looks much better. 
However, at an even higher Mach number (Fig. 11c), 
the calculations predict values considerably below the 
measured results for low Reynolds number (viscous 
effects greater than inviscid), while for high Reynolds 
number, the results are satisfactory. For the highest 
sweep angle tested, 70° (Fig. 12), the calculations give a 
pressure distribution slope which is considerably higher 
than experimentally measured and, in general, predict 
pressures much too high over the wing, particularly up- 
stream of the s/t = 40 to 60 station. A similar trend is 
shown at higher Mach numbers. 

A summary of the effect of varying sweep angle with 
constant Ke, , is shown in Fig. 13 for Mach numbers of 7 


and 17. Ata Mach number of 7 the pressure distribu- 
tion slopes vary from about —0.64 to —0.56 with a 
major divergence to —0.23 for A = 70°. For a Mach 


number of 17, the slope of the curves for A = 0° to A = 
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60° vary from —0.56 to —0.52, with the value of —0.42 
at A = 70°. 

It is quite noticeable that for a sweep angle of 70°, the 
pressure distributions have a different variation than 
for the other sweep angles. There is a rather sharp 
break between 60° and 70°. The off-axis static pres- 
sure orifices were carefully examined to see whether 
there was any indication of tip effects. However, the 
spanwise results for the 70° swept wing did not deviate 
significantly from the results obtained at the other 
sweep angles. There is still, however, some possibility 
that the apex effects are important for this test; but 
further studies will be needed to verify this suspicion. 

One further piece of interesting information can be 
obtained by examining the pressure at a given station 
(for a given Mach number and leading-edge thickness) 
as the sweep angle is varied. A typical set of results for 
several s/t’s at a Mach number 17 is shown for the 
highest Reynolds number tested in Fig. 14. The slope 
of the straight lines drawn through the points is about 
1.1 as compared to the geometrically predicted 4/3 
power. There is some definite indication that, for the 
highest sweep angle tested, the pressures are somewhat 
higher than predicted by this linear power variation. 


Concluding Remarks 


Before summarizing the results which have been pre- 
sented, several important points should be made with 
regards to the calculations and the tests presented 
herein. The calculation technique of simply adding the 
viscous and inviscid terms cannot be correct for most of 
the cases reported herein where the viscous and inviscid 
terms are about the same order of magnitude. The 
idea of adding two such effects, without any interaction, 
can only be considered a very crude picture or perhaps a 
way to simply describe the experiments. Probably the 
major restriction on the present study has been the 
implied two-dimensionality of the flow. Such two- 





SWEEP ANGLE 979 


dimensional flows were easily obtained for the unswept 
wing; but the swept wing requires a closer examina- 
tion. One criterion of two-dimensionality is that the 
pressure is constant along lines parallel to the leading 
edge. This was borne out by the pressure orifices 
located off the main line of orifices. However, this still 
does not definitely prove that results are free of any 
apex effect. A rather complete study of this apex effect 
is currently being carried out; and at present, it is 
believed that the results presented are free of such end 
effects. 

The actual flow over a swept wing is definitely not of 
such a two-dimensional character, and this probably can 
be best shown by an oil streak photograph (Fig. 15). 
At least on the surface, the flow direction is inclined at 
an angle of about 15° to the free stream and, at other 
Mach numbers and sweep angles, flow direction angles 
Although the 
“pressure constant along lines parallel to the leading 


as high as 25° have been measured. 


edge’’ satisfies the simple two-dimensional criterion, 
this pressure field generates a lateral pressure gradient 
which influences the boundary layer. This, in turn, 
can have a considerable effect on the developrrent of the 
pressure distribution. Considerable further work is 
being carried out on these three-dimensional flows. 

With the above observations as a background, the 
results of the present studies, carried out in the Mach 
number range of 7 to 18 for leading-edge Reynolds 
numbers in the range of 3,000 to 25,000 and for s/t 
dimensions up to about 200, can be summarized as 
follows: 

(1) For the unswept wing, the simple addition of the 
calculated viscous and inviscid effects agrees with the 
measured pressure distributions over the entire range 


tested to within about 15 percent. The viscous correla- 








Fic. 15. Oil streak photographs on a swept plate, ¢ ~ 0.03 inches, 
A = 45°. Top: M = 8, bottom: M = 17. 








980 JOURNAL OF THE AEROSPACE SCIENCES—DECEMBER 1961 


tion for variable leading-edge thickness was within the 
order of +5 percent. The inviscid Mach number 
dependence, as inferred by subtracting the theoretical 
viscous term, was found to vary as the 2.2 power as 
compared to the theoretical 2.0 power. Using the 
results herein and other available data in the Mach 
number range of 17 to 24, for a fixed Re, , of about 
7,000, the Mach number dependence for the combined 
viscous-inviscid effect was found to be about 2.4. 

(2) For the swept wings, it was found that the simple 
addition of the calculated viscous pressure distribution 
(based on the free-stream conditions with no effect of 
sweep), and the inviscid pressure distribution (based 
on blast wave analogy and the normal Mach number) 
does not give consistently satisfactory results. Devia- 
tions above and below the calculated values by as much 
as 40-50 percent were measured over a range where the 
viscous effects varied from about 40 to 140 percent of 
the inviscid pressures calculated. 

(3) The results for the high sweep angle of 70° 
appear to deviate rather sharply from the usual results 
at lower sweep angles. The pressure distributions were 
considerably different than predicted, particularly in 
slope, with some of the measured results as much as 50 
percent less than the calculated values. 

(4) The effect of sweep at a constant Reynolds 
number gives a pressure distribution which varies with 
distance along the wing with an exponent between 
about —0.53 and —0.58, with a rather strong decrease 
occurring for the 70° swept case. 

(5) The pressure at a given station, fora fixed Mach 
number and a given leading-edge thickness, varies as the 
cosine of the sweep angle to the 1.1 power is compared 
to the 1.3 power predicted from general geometrical 
sweep considerations. 

It is obvious from the presented results that the 
simple pressure distribution over a swept blunt two- 
dimensional wing at zero angle of attack cannot be 
calculated with sufficient accuracy to carry out studies 
of hypersonic wings. Considerable further theoretical 
work is required, as well as additional experimental 
studies, to determine the main parameters of this flow. 
Three-dimensional effects will obviously have to be 
investigated and will be a major contribution to the 
actual flow over delta wings. 
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Nonuniformities in Magnetogasdynamic 
Channel Flow 


Edward T. Pitkin* 


The Marquardt Corporation, 
Van Nuys, Calif 


March 6, 1961 


“s QUASI-ONE-DIMENSIONAL tmagnetogasdynamic channel 
flow equations have been presented! for the case of mutually 
fields As it 


pressure, velocity, 


perpendicular applied electric and magnetic was 
later pointed out,? the flow variables—e.g., 
current, field, body 
to be taken as average values, 
tions arises only in evaluating the induced part of the magnetic 
field by Biot-Savart This induced field is 
then subtracted total 


ascertain the 


electromagnetic force, etc.—are 


and the need for Maxwell's equa 


magnetic 


means of the law 


field used in the channel-flow 
field to be 


calculation, the 


from the 


solution, in order to applied at each 


position. Following the indicated resulting 


expressions for the x and zs components of the induced field are 
seen to exhibit a dependence not only on position, but on channel 
height, a, width, 6, and length, Z. The notation is that of 


Ref. 2 with x indicating the flow direction 


ad ! 
he = f j,(é) In X 
- a 
iiy+0+ Vie £)2 + (y + 5)? + (2 a)? 
« = = x 
Uy —b + V(x — £2 + (y — 0)? 4 (zs — a)?] 
[ly —b + V (x — £)2? + (y — 5)? + (2 + @)? 
dé (1 
ly + h+ WV (x ~£)2 + (y +b)? + (s +a) 
PL 
h. | jy(é) X 
0 
J (s + a)(y b 
tan™! a 
' yr —i)V(x —f)? +(y —b? + +4 
(2 — a)ly 5) 
tan™ + 
(x — £)V (x — &)? + (y — 5)? + (2 — a)? 
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: (zs + a)(y + d) 
tan r _ 
(x — —)V (x — &)2 + (y +b)? + (2 + a)? 
(g —a)(y + db) 
tan™! i : t dé (2) 


(x — V(x — &)2? + (y + bd)? + (2 - a)?4 


By again taking an average value of h(x) for given channel 
dimensions one can however obtain a useful approximation to 
the ideal field distribution. This is certainly in order since the 
total field, H.(x), has already been taken as an average 

The x component of the induced field will generally be difficult 
to nullify by an applied field, so h, H,, the total field. The 
body forces j7,H, and the nonuniform 7,/7, will then cause a 
distortion of the velocity profile. It may be shown that the 
j,H, force predominates in long slender channels.?__ The resulting 
distortion in the velocity profile may then be estimated. Con- 
sider a slender channel in which £,(x) and H,(x) vary slowly so 
that the results of Ref. 2 apply—.e., 


OH,/0z > OH,/Ox (3) 
Consider only variations in the z direction to estimate the change 
in the velocity: 
OH,/0s = 4n0(E, — uz) (4) 
and the z component of the momentum equation with w < u 
gives 
Op/dz = o(E, — uH.)H, = (H,/41r)(0H,/0z) (5 


The velocity variation referenced to the centerline value is 


bu 1 ff; Ow... : 
= - & (6) 
ut. ut J9 Of 


If the flow is initially uniform, the Bernoulli law will give a good 
approximation to the pressure-velocity relationship in the trans- 
verse direction: 

‘ 


ou/Op = —1/pu ( 
Letting pw remain constant at the centerline value allows an 
iterative integration of Eq. (6) with the aid of Eqs. (5) and (7 
6u/ue. = |H.2/87r(pu?)¢] (H,/H-)? (8 
A similar integration of Eq. (4) with uw held constant for the first 
approximation gives 
H,/H, = 4routs|(E,/ujH.z) — 1] (9 
Combine Eggs. (8) and (9) and define three nondimensional param 
eters: 4rou¢z 4 R,,, the magnetic Reynolds number based 
upon z; (uv 4 wp)¢./H,AA, the Alfvén number; and u¢.H,/E,4n,, 
the instantaneous conversion efficiency. Then 
bu/ue, = 1/2} [(1 — 9.)/ne] (Rm/A) $? (10) 
Each of the nondimensional parameters is based upon centerline 
average values and can be estimated readily. As long as the 
right-hand side of Eq. (10) < 1, nonuniformities in the velocity 


profile may be neglected. 
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+ 
Comment on ‘‘On Column Behavior’’ 


lain Finnie* 
Shell Development Company, Emeryville, Calif. 
April 17, 1961 


T° A RECENT ARTICLE! in the Reader’s Forum, S. Sergev dis- 
cusses the behavior of pressurized slender columns with 


* Presently at the University of California, Berkeley 4, Calif. 
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either open ends or closed ends. It should be pointed out that 
the solution to this problem was given by J. A. Haringx? in 1951, 
The original article was in Dutch but was translated into English 
in The Engineer's Digest? and was reviewed under the heading 
“buckling problems” in Applied Mechanics Reviews.* 

In his article Haringx makes the interesting observation that 
the general buckling equations in the texts of Biezeno-Grammel, 
Love, and Timoshenko (elastic stability) do not predict the 
instability of the pressurized tube with open ends. The texts of 
Fliigge, Pfliiger, and Prescott on the other hand contain general 
equations which do predict this type of instability. As Haringx 
points out, the difference between Fliigge and Biezeno-Grammel 
is that the latter authors did not include the axial Component of 
pressure in the calculation of the additional load due to hydro 
static pressure in the deflected condition 

It may be of interest to note that this type of buckling is recog- 
nized in the design of oil-well casing. This was discussed by 
Woods? in 1951, and illustrated by a model. Since that time this 
particular form of instability has been quite well-known in the 


oil industry. 
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Internal Pressure, The Engineer's Digest, Vol. 12, No 
Oct. 1951 

‘ Biezeno, C. B., Review of 
to Internal Pressure,’’ by J. A. Haringx, Applied Mechanics Reviews, Vol. 4, 
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‘Instability of Thin-Walled Cylinders Subjected 


Numerical and Approximate Solutions of Re- 
Entry at Large Angles of Inclination 


W. H. T. Loh* : 
Vought Research Center, Chance Vought Corporation, 


Dallas, Texas 
February 20, 1961 


— ggomaed INTEREST in orbital re-entry into a planetary at 
mosphere has resulted in some approximate analytical solu- 
tions. For ballistic-type re-entry without lift at sufficiently large 
angles of inclination where both the gravity force and centrif- 
ugal force are being neglected, the approximate solutions of 
Gazley,! Allen and Eggers,? and Chapman? are available. For 
gliding-type re-entry at positive lift/drag ratio and small angles 
of inclination, the approximate solutions of Allen and Eggers‘ 
are available. For re-entry of a satellite with a small lift/drag 
ratio or re-entry of an orbiting vehicle with a small initial angle 
of inclination, the approximate solutions of Chapman? are avail- 
able. For tactical and strategic reasons or for avoiding natural 
radiation hazards, quick descent at large angles of inclination is 
also being considered. For this type of re-entry approximate 
solutions at constant negative lift/drag ratio and large angles of 
inclination were given previously.° Since no other analytical 
or numerical solutions were found in the literature at that time, 
a case of nearly zero lift/drag ratio was used for order-of-magni- 
tude comparisons. However, it was pointed out there® that 
the approximate solution has a singular point at zero lift/drag 
ratio and therefore it is to be expected that the solution probably 
will yield a better result at negative lift/drag ratio greater than 
0.1. Recently, exact numerical machine solutions for re-entry 
at constant negative lift/drag ratio and large angles of inclination 
were obtained at the Mathematics Institute of the Federal 
Technical Institute of Zurich, while the author was there. A 
comparison with the approximate solutions given previously® 


* Staff Assistant to the Associate Director of the Research Center. 
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Fic. 1. Comparison of analytical approximate solution with exact numerical machine calculations 


shows good agreements for negative lift/drag ratio greater than 
).1. Since both the numerical solutions and the comparisons were 


not known previously, it would seem to be of interest to present 
1 


them both here (Fig. ! 
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On Computing Inlet Entropy Increases for 
Nondissociating Gases 


C. W. Hartsell 
The Aerospace Corporation, Los Angeles, Calif. 


January 31, 1961 


SYMBOLS 


Cy = specific heat at constant pressure 
VW = Mach number 

P = pressure 

R gas constant 

Ss entropy 

7 = temperature 

p = density 


¥ = ratio of specific heats 
6 constant, for air @ = 5,500°R 


Subscripts 
free-stream conditions 
= conditions at the end of compression 
} = perfect gas 
1 = total or stagnation conditions 


A“ ANALYSIS of gasdynamic problems in the intermediate 
region between perfect-gas techniques and real-gas calcu- 
lations using Mollier diagrams has been developed.! The ratio 


of specific heats is allowed to vary as a function of temperature 
along with the assumption of zero dissociation. It is the pur 
pose of this note to present a simple means for using the tech 
nique of Ref. 1 in inlet-performance calculations. The parameter 
of basic interest in inlet calculations is the entropy rise during 
the compression process. If the entropy rise is known the inlet 
efficiency may be easily evaluated. Thus an equation presenting 
the entropy rise as a function of conditions at the end of com 
pression plus the free-stream conditions would be a great help 
to the inlet analyst 

Combining the first and second laws of thermodynamics for a 


thermally perfect gas—i.e., P pRT-—we have 


fas J C,(dT/7 Rf dP/P (I 
1/R fas f ly/(y — 1)])(d7/T { ap gs 2 


For a nondissociating gas may be related to temperature as 





follows: 
1 ae 3 
. 4 3) 
j 1 + (yp — 1) {(9/T)? le/7 /(e8/7 — 1)3]} 
Insert Eq. (3) into Eq. (2) for y, 7/5 
1 7 (aT (*)' o/T dT “dP , 
S = oa oes (4) 
ae “453 ie i +See P 
integrate between stations (@) and (2 
AS 7 ee est 
: In 72/7. — In P2/P. + In [(e 
R 2 
0 6 6/T>» 6/T 
1)/(e9/72 — 1)] + — + (5 
j T. T x) T: 1 6/7 1 
For use where 7... is low, such as at free-stream conditions—.e., 
e9/T= > 10%—we may approximate Eq. (5) as follows: 
AS 7 Le 
= - In 72/7 -In P2/P. - 
R 2 
0 O/T 
In (e9/7? — 1) 4  & 6/T2 l " 
2 \e = 


A comparison between the entropy increases predicted by 
Eq. (6) and the results from perfect-gas and real-gas solutions is 
presented in Fig. 1 for the case of a normal shock. It may be 


seen that the present solution satisfactorily bridges the gap be- 
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Fic. 1. Normal shock entropy increase. 


tween the two techniques. For this particular case it is seen 
that the usefulness of the present technique is approximately 
4 < M.. < 8 which adequately covers the present area of interest 
for inlet designers. It must be mentioned that the choice of a 
normal-shock test case is a very severe test, as higher efficiency 
inlets delay the onset of dissociation extending the range of use- 
fulness of this technique. In this respect the technique of Ref. 1 
provides very good values of Py.. and T,7., up to hypersonic 
speeds at all reasonable altitudes because of the negligible dissoci- 


ation present. 
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Further Comments on ‘‘A Drag Hypothesis for 
Jet-Flapped Wings’’ 


G. K. Korbacher 

Associate Professor of Aeronautical Engineering 
Institute of Aerophysics, University of Toronto 
Toronto, Ontario, Canada 

June 22, 1961 


ie HIS AUTHOR’S REPLY, Professor Reid! took exception to a 
statement in my recent Readers’ Forum note! which says 
that Foley’s and Reid’s ‘‘finding is substantiated by test results 
of NGTE? * and ONERA.’’* Reid claims that ‘‘Foley’s finding 
has been substantiated only by the results of subsequent——and 
not yet published—Stanford work.’’ In this writer’s opinion, 
however, his claim is not substantiated by the facts. 

The facts are that when the many-years-older experimental 
data of NGTE (Dimmock? *) and ONERA* are plotted in the 
same way as those of Foley and Reid® (their Fig. 3), they not 
only substantiate Foley’s and Reid’s finding of 94° 7 thrust re- 
covery, but even better it. This is shown in Figs. 1 and 2 (from 
Ref. 6), which are a replot of NGTE’s and ONERA’s test results. 

Foley’s and Reid’s contribution is thus seen to lie in the novel 
way of analyzing and plotting the experimental data (see Fig. 
3 of Ref. 5) such that the thrust recovery is presented as the slope 
of curves at constant lift coefficient. Their experiments, how- 
ever, are in no way unique. They support the earlier English 
and French work but definitely do not supersede it. 
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Fic. 1. Drag variation at constant lift (Dimmock’s results for 
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Author’s Final Reply 


Elliott G. Reid 
Dept. of Aeronautical Engineering 
Stanford University, Stanford, Calif. 


July 31, 1961 


r i SHE WRITER must categorically reject Professor Korbacher’s 
opinion on the basis of erroneous interpretation of data, 
unanimously contrary opinions of the cited British work by those 
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who did it, and concurrence in the appraisal of GNERA results 
as nondefinitive 

The error of interpretation is that of accepting the slopes of 
constant-C;, curves of Cp vs. Cy as the magnitudes of thrust- 
recovery fractions. To be sure, ideally complete thrust recovery 
at fixed Cy, would necessarily be characterized by dCp/dCy = 
-1.0 but satisfaction of this condition, alone, is not sufficient to 
prove the achievement of such recovery. Convincing evidence 
of this fact is provided by the following results of tests of a modi 
fied version of Foley’s model: At Cy, = 2.0, dCp/dCy ~ —1.7, 
but the thrust recovery at Cy = 0.5 was less than 65 percent. It 
must be emphasized that the only irrefutable proof of thrust 
ecovery is drag reduction. 

With no little chagrin, the writer must confess to at least par- 
tial responsibility for Professor Korbacher’s confusion which un- 
doubtedly stems from the carelessly erroneous statement in the 
1959 Forum article® that, ‘‘Thus, Fig. 3 clearly indicates that a 
large and substantially uniform percentage of the calculated jet 
reaction was transformed into upstream thrust etc.’’ Actually, 
the stated fact was substantiated only by Fig. 4—since a fifth 
(more readily interpreted) chart was then deleted for want of 
space. That figure is reproduced below in the interest of clari- 


fication. 
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Professor Korbacher’s contention that complete thrust re- 
covery was demonstrated by early NGTE work constitutes a 
lone dissent from all qualified British opinion known to this 
writer. Dimmock’s contrary appraisal of his own work was 
cited in the Author’s Reply of May 196i. In the most recent 
summary of British jet flap work,! Williams, Butler, and Wood 
state that ‘‘the sectional thrust lies between the corresponding 
horizontal component and full jet momentum.”’ 

Even more puzzling is the attempt to discredit the originality 
of Foley’s finding by the citation of earlier ONERA work. An 
addendum to the legend beneath Fig. 5 of the critic’s own re- 
port® reads as follows: ‘(QNERA Test Results—Due to Lack 
of Information—Cannot be Converted to a Cpr vs. Crzr? Plot, 
A Fact which Makes this Figure Useless for Quantitative Com- 
parisons.)’’ As was stated in the Reply of May 1961, Foley’s 
investigation was largely motivated by concurrence in this 


appraisal. 
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An Aeroelastic Analysis Procedure for 
Flap-Type Control Surfaces} 


William S. Johnson, Jr.* 
General Dynamics/Fort Worth 
March 21, 1961 


SUMMARY 
This note describes a steady-state aeroelastic solution amenable 
to all types of planform configurations and aerodynamic, struc 
tural representations. The solution considers aspects of the 
problem of limited coverage in the literature.! 


SYMBOLS 


CHs hinge-moment coeflicient due to aileron deflection 
Cis = rolling-moment coefficient due to aileron deflection 
( ly rolling-moment coefficient due to damping in roll 
D = plate flexural stiffness 
qd = free-stream dynamic pressure 
, Clg [{“% 
¢ = rolling effectiveness 
( lp I ( ly R 
Subscripts 
I = elastic 
I flap surface 
R = rigid 
Mu main surface 


a SOLUTION for aeroelastic equilibrium is possible on 
the digital computer with the techniques of matrix algebra 
Reviewing schematically: 

From a structural subprogram a stiffness-matrix equation is 
obtained: 


[stiffness matrix] X } deflections} }loads} 


From an aerodynamic subprogram, an aerodynamic influence 
matrix: 


; ; ‘ ; ; 
q |aerodynamic matrix] X {deflections} }induced loads} 


For aeroelastic equilibrium: 
[stiffness matrix] X } equilibrium deflections} ) }rigid loads} + 


i 
induced loads} | 


or 


[stiffness matrix] — g [aerodynamic matrix]] 


}equilibrium deflections} ) rigid loads} 


and 


}equilibrium deflections} [stiffness matrix 


g |aerodynamic matrix]]~!  }rigid loads} 


The individual arrays of this equation may be constructed to 
include all aspects of the problem: main and flap flexibility, 
structural and aerodynamic interaction of the main and flap 
surface, and the flap-actuator system. The general form of the 
aeroelastic matrix equation with provisions for all of the above 


variables is: 


[B}"} [[A) [Su] (M0) + [DI)7{HL)7 [Se] [HL] [DI)] 


l O 
g[( A)": (DI) "|HL}? [G] [ZS] [7C| [A] 77 B| x 
0 DI 
\e/ \ \L! ; ‘ke & 
. , = [BB]? M\" +. + [DI|\T[|HL|Ts .. 
la\mu [B| ( (7\ (T\ r\ 
pas p’ 


Each variable is represented as an independent matrix array and 
may be varied arbitrarily for parameter studies. The logic of 


t+ This paper is based on a thesis of the same title submitted to Southern 
Methodist University in 1959 as partial fulfillment of the requirement for 
the degree of Master of Science in Mechanical Engineering. Thesis advisor 
was Dr. J. S. Archer, now of Space Technology Laboratories. Inc 


* Project Structures Engineer. 
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the equation can best be understood by thinking of it as a simple 
algebraic equation in which the matrix arrays are coefficients 
performing a specific function. 


STIFFNESS MATRIX 


[Su], [Sr] 

The main- and flap-surface stiffness matrices are obtained 
independently of a boundary fixity by any of the available 
methods. Orientation is generally to a deflection and/or 
rotation (6, a) displacement, and a shear and/or torque (L, T) 





load: 
. fa) (LZ) = 
Sad) ob ae ~ (Tha (2) 
fe) (L) ? 
sar} = {=} 
ale (T\r - 
[iM] 


This matrix contracts the displacement column of the main 
surface-flap combination to that of the main surface: 


{a} {o) 

a. «fear § ( 

la\m 4) 4) 
{a} p’ 


[DI] 

The flap surface is flexible rotationally with respect to the 
main surface. If the flap is actuated by one or more position- 
sensitive rams, the flap is ‘‘fixed’’ rotationally with respect to 
the main surface at the ram stations. If the flap is actuated by 
two or more hinge-moment-sensitive rams, a definition of aileron 
deflection for the elastic structure must be established. In 
such a system, control-surface deflection is usually read at one 
inboard flap station. This definition of flap position is projected 
into the aeroelastic equation as an equivalent rotational ‘‘fixity”’ 
of the main and flap surfaces at the deflection-indicator station: 


arp = ay (5) 
For all arrangements: 
\@ 0 
a an ri} l (6) 
la\w a\ au 
falter { a} e’ 


|HL| 

If the flap surface is of a size and stiffness to contribute bending 
restraint to the main surface, this structural interaction can be 
considered by stipulating that vertical deflections of the two 
surfaces at the hinge-line be equal: 


{o) (oe) oes 
el in AL} 4) 
lalr 


AERODYNAMIC MATRIX 


[A] 
The aerodynamic influence matrix is of such a form that in- 
duced pressures at prescribed points can be calculated knowing 


local slopes or mode coefficients. 
[IT] 


This interpolation matrix calculates slopes at the aerodynamic- 
matrix pattern from deflections and/or slopes at the stiffness- 
matrix patterns. If the aerodynamic matrix is constructed for 
various twist and/or camber modes, the [77] matrix may be a 
least-squares or similar operator for calculating the required 
coefficients. 

The pressures must be integrated over suitable areas to obtain 
local shears and torques compatible to the stiffness matrices 
load points. For convenience, this operation is accomplished 


in two steps: 


. UC]  chordwise 
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Fic. 1. Influence of the control system on aileron effective- 
ness. Semiplan exposed aspect ratio = 1.55, Mach number = 
2.0. (Left) Aileron rolling moment per unit hinge moment. 
(Right) Rolling effectiveness. 


and 


[ZS] spanwise 


[G] 

This matrix is peculiar to a system with two or more hinge- 
moment-sensitive rams. It imposes the condition that induced 
torsional loads on the flap must load the main surface and react 
the flap surface according to the sizing of the actuator rams. 
Wing thickness usually limits actuator size: the inboard rams 
carry a larger portion of the hinge-moment than the outboard 


rams. 
EQUILIBRIUM DEFLECTIONS 


[B] 

The boundary-condition matrix imposes a symmetric, anti- 
symmetric, or arbitrary support condition on the main-surface 
stiffness matrix: 


\e] \o/ 
19} ue = [B?t (8 
la) la\a 
‘=. ee 
LQsPr LQsF 


Ricip Loaps 
\L( \LI 
iT\a (T\r 
The rigid shear and torque loads are calculated for the com- 
plete stiffness node pattern, Eqs. (2) and (3). 


SAMPLE PROBLEM 


To illustrate the significance of the actuator-system charac- 
teristics in a control effectiveness study, a simple, constant- 
thickness, plate wing was analyzed for various ram arrangements: 


A = a “plain” surface analysis in which the flap is assumed 
to twist as a simple extension of the main surface 

B = single actuator at the inboard end of the flap 

C = single actuator at the outboard end of the flap 

D = two actuators, one at each end of the flap, each reacting 


half of the hinge-moment 


Main- and flap-surface stiffness matrices were calculated by 
the linear solution of Stein-Sanders.2 A full-fixity boundary 
The flap- 

Results 


condition was imposed at the fuselage junction 
bending restraint was considered in all arrangements 
are shown in Fig. 1. 


CONCLUSION 
No direct solution for an optimum control system (actuator 
sizing and arrangement, main/flap stiffness) is evident. An 
acceptable arrangement can be determined only from a systematic 
study of these variables across the flight envelope. Results are 
frequently unpredictable: stiffening a flap can result in less 
control effectiveness; the best configuration subsonically may 


not be the best supersonically. The suggested type of analysis 
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is easily covered with the matrix solution of Eq. (1). 
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A Simplified Approximate Method for the 
Calculation of the Pressure Around Conical 
Bodies of Arbitrary Shape in Supersonic and 
Hypersonic Flow 


Willi F. Jacobs 
Lockheed Aircraft Corporation, Georgia Division 
March 3, 1961 


te CONICAL-FLOW SOLUTIONS are available only for circular 
cones at zero angle of attack.! For nonaxisymmetric cones 
or cones at angle of attack, only approximate methods exist. 
These methods are generally quite complicated and further 
limited to certain body shapes or certain Mach-number ranges.?, 4 
A great need was therefore felt for a simple approximate method 
applicable to any arbitrarily shaped conical body at zero inci- 
dence as well as at angle of attack. 

Such a method has been developed recently at Lockheed® and 
is presented here in abbreviated form. The method is based on 
the ‘‘equivalent-cone’”’ theory.‘ This theory determines the 
pressure on a conical body utilizing information for a symmetric 
cone at zero angle of attack with the same normal component of 
the free stream with respect to the surface as the local element 
of the body considered. This method works relatively well at 
high Mach numbers. However, it is quite inconsistent at lower 
Mach numbers, especially for bodies which deviate considerably 
from circular cones. The equivalent-cone method does not give 
satisfactory results, mainly due to the fact that it considers only 
the local surface element on the body independent of the other 
body elements in the Newtonian-theory manner. 

These interference effects between the different surface ele- 
ments have been taken into account in the newly developed 
method by adding to the equivalent-cone value of pressure an 
additional correction term for the pressure exchange along the 
surface of the body. This new method may be called the ‘‘im- 
proved equivalent-cone method.’”’ It covers the whole super- 
sonic Mach number range for which the conical flow conditions 
with attached shock at the apex of the body are fulfilled. 

The assumptions made for the correction term are that the 
pressure exchange along the surface is proportional to the pressure 
difference between the equivalent cone and circular cone. This 
reference circular cone is a cone having a pressure equal to the 











Fic. 1. Symbols and body geometry. 
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average pressure of the equivalent body. This correction term 


is further dependent on the Mach number in the conical flow 
field. The pressure coefficient on the conical body of arbitrary 
cross section is then given by the expression: 


Cp = (P — po)/g = Cre — (Coe — Coe) f(M 
where C,, is the pressure coefficient on the equivalent cone and 
Cpe its average value over the body. It has further been assumed 
that the Mach-number correction, f(.1/), is of the form: 


f(M) = 1/M" 


This satisfies the flow conditions both at very high Mach numbers 


where the pressure is given by the value of the equivalent cone: 
(Cr)m.1 = Cp 


and at low Mach numbers for which the pressure approaches the 
constant value of the circular cone: 
(Cp)mM—+1 = Cp 

The exponent » in the Mach-number correction has still to 
be determined. The assumption that the pressure exchange 
along the surface of the body is proportional to the Mach angle 
in the flow field behind the shock wave leads approximately to 
n= 1,oFr: 


f(M) = 1/M" = 1/M, 


where J/, is the Mach number on the circular reference cone 
with the average pressure of the equivalent cone. The final 
formula for the ‘‘improved equivalent-cone method”’ is then: 


Cp = Coe — (Coe — Coe)(1/Me) = Cell — (ACp/Cpe)(1/Me)] 


A comparison of the new method with available experiments 
on elliptical cones and cones at angle of attack gives very good 
agreement. This is illustrated with a few examples. Fig. 1 
shows the symbols and body geometry. The pressure coefficient, 
Cp, is plotted against the peripheral angle @ of the body cross 
section. Fig. 2 presents the results for an elliptical cone with 
the ratio of minor to major axis e = b/a 1/3 and the slender- 
ness ¢ = x/V S; = 2.93 at a Mach number of 1.81.2. The com 
parison shows that the ‘“‘equivalent-cone’’ theory results in pres- 
sures in the high-pressure region which are too high and pressures 
in the low-pressure region which are excessively low. However, 
the new method is in complete agreement with the measurements 
Similar results are shown in Fig. 2 for an elliptical cone with e = 
0.56 and o = 1.875 at a Mach number of 6.04.4 The applica- 
bility of the new method for conical bodies at angle of attack 
is proved in Fig. 3 for a circular cone with semi-apex angle @, = 
11.42° at @ 5° and at a Mach number M = 5.05.4 

Summarizing, it can be said that very good agreement has 
been obtained between results of the new method and available 


measurements, 
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Comments on ‘‘On the Calculation of Natural 
Modes of Free-Free Structures’”’ 


David Williams 
Deputy Chief Scientific Officer, Structures Dept., 
Royal Aircraft Establishment, Farnborough, Hants., England 


March 6, 1961 


r HIS NOTE! John Dugundji sets out to ‘‘indicate certain 
interesting relations which may not have been generally ap- 
parent heretofore,”’ in the belief that the points he makes have 
not previously been made. May I point out however that in 
this he is not correct. In my recent book? the points he raises 
are discussed in detail—S4.14, 4.15 consider the case of the 
simple free-free beam and the whole of Chapter 12 is devoted 
to the three-dimensional problem of deriving the normal modes 
of vibration of a complete airplane structure. 
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Author’s Reply 


John Dugundji 
Assistant Professor of Aeronautics and Astronautics, 
Massachusetts Institute of Technology, Cambridge 39, Mass. 


April 14, 1961 


M* READERS’ FORUM NOTE! dealt specifically with some special 
features of the flexibility influence coefficient matrix 
method of calculating the natural modes of free-free structures. 
In this method, attention is called to the complete independence 
of the [R] matrix from the particular set of restrained-structure 
flexibility influence coefficient matrix [C] chosen to describe the 
structure. This then resulted in the interesting phenomenon 
that totally different restrained-structure dynamic matrices 
[C] [m], possessing different eigenvalues and eigenfunctions, when 
premultiplied by a numerically identical matrix [R], would then 
all possess the same eigenvalues and eigenfunctions (the free- 
free frequencies and modes). 

There are several other methods of calculating natural modes 
of free-free structures. Dr. Williams in his book? treats this 
problem in considerable detail with a finite difference approach. 
The above points regarding the [C] and [R] matrices do not come 
up explicitly there. 
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A Semiempirical Description of the Discharge 
Characteristics of the Converging Section of a 
Low-Density Hypersonic Nozzle* 


Frederick O. Smetana 


Research Scientist, Engineering Center, 
University of Southern California, Los Angeles 7, Calif. 


March 10, 1961 


—. DEVELOPMENT of low-density hypersonic wind tunnels 
has created a number of new nozzle design problems. It 
has been observed, for example, at USCEC that when the up 
stream stagnation pressure of a convergent-divergent nozzle is 
decreased below a certain value the flow at the nozzle exit is sub- 
sonic no matter how low the exhaust pressure. One possible ex- 
planation for this is that the flow at the throat never becomes sonic 
because of large viscous effects. As shown below, this effect 
can be important if the stagnation pressure is low or the throat 
diameter is small. Both of these conditions are usually present 
in low-density wind-tunnel nozzles because of the low mass flows 
which can be accommodated by the pumping systems 

Consider the convergent portion of an axisymmetric nozzle as 
shown in Sketch A. When 102 < Reg < 10‘, one is probably 





SKETCH A, 


justified in assuming that the flow at the throat consists of an 
isentropic core and a laminar boundary layer with displacement 
thickness 6*. The discharge coefficient, Cz, can be represented, 
therefore, by 


Cu = [(d — 26*)/d]? ~ 1/[1 + (46*/d)] (1) 


If one assumes that 6* = (K,s)/V Re,, and s/d = constant, Ko, 


one may write 
6*/d = (KiV K2)/V Rea 


Hence, Eq. (1) can be written 

Ca = 1/{1 + [(4KiV K2)/V Real} = 1/[1 + (B/W Rea] (2) 
The value of A, depends upon dP/dx. One may estimate the 
value of AK, from the exact analysis of flow over a cylinder. In 


that analysis 
6* = 0.8R/V Ua2R/v 


at the 90° point. For R/d = 1.66 ~ Ko, B = 4.12. Tests ofa 
nozzle with this entrance geometry gave B = 4.0. The calibra- 
tion data are shown in Fig. 1. Dry nitrogen was the test fluid. 

Eq. (2) can be expected to yield reasonable agreement with 
experiment down to Reynolds numbers for which the boundary 
layer completely fills the throat. If one assumes that when this 
occurs, 6*/6 = 0.27, then one has Cy = 0.65 which agrees with 
experiment. 

As the Reynolds number is decreased below this point, the 
flow velocity at the throat decreases. This comes about be- 
cause the pressure at the throat does not decrease—as near as can 
be determined experimentally—compared with the upstream 
pressure while the viscous dissipation continues to increase. 
For such flows, the discharge is the same as for incompressible, 
viscous flow for which m ~ AP in a straight pipe. Isentropic 
nozzle theory, however, gives m ~ V AP for constant P. so that 





* The research reported in this paper has been sponsored by the Wright 
Air Development Division, Air Research and Development Command, 


United States Air Force under Contract AF 33(616)-7102 
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Fic. 1. Discharge characteristics of Venturi meter. (1) 


Three different McLeod gauges were used. The divergence of 

data at the left of the graph is due to remaining differences in 

indications of two types of McLeod gauges. (2) The curves 
f 


shown are represented by the equations C, {1/1 + (4V Re)] 
Ca > 0.65; Ca = (VRe/11.4) Ca < 0.65 


in a slowly contracting channel 7/Minviseia ~ WV Re Thus, 
Ca = V Rea/C (3 
From the viscous pipe-flow equation and the isentropic nozzle 
equation, one estimates a value for C = 7.14 for nitrogen. In 
a rapidly converging channel there is more dissipation and hence 
the value of C should be greater. The data in Fig. 1 indicate 
that 11.4 is a good value for C. This is 2.86 times B for the same 
nozzle. <A check of two sets of data in the literature on incom- 
pressible values for Ca vs. Reg indicates that this ratio may be 
quite universal. If this is the case, then one may describe the 
discharge characteristics of a given nozzle in terms of a single 
constant which depends on geometry: 
, l c , ‘- ‘ 
Ci = for Ca => 0.65; Ca = 
1+ (B/V Rea 
V Reg. ’ a 
> for Cy 4 0.65 (4 
2.86 B 
On the basis of this interpretation, some rather interesting 
results may be obtained pertinent to the design of nozzles for low- 
density wind tunnels: (1) When B = 4, the minimum stagna- 
tion pressure for which the flow will accelerate isentropically to 
supersonic velocities downstream of the throat is 
Pmin © (0.688/d)(300/T,)'-* (5) 


where P is in mm. Hg, d in cm, T, in °K, and the gas is No; (2) 
nozzles with very gradual entrances will have high values of 
s/d and therefore of B which leads to a higher P,,,,. Thus, the 
nozzle entrance should be the sharpest possible which avoids 


separation. 


One-Dimensional Heat Conduction With 
Arbitrary Heating Rate and Variable Properties 


J.C. Y. Koh 

Research Specialist, Aero Space Division, Boeing Company, 
Seattle, Wash. 

May 1, 1961 


SYMBOLS 
= specific heat 
= thermal conductivity 
thickness of slab 
= surface heat flux 
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‘ 
— =x Fic. 1. Physical model. 
. "i ‘ 
ms 
t = time 
t& = penetration time 
Ue surface temperature from exact solution 
t = Eq. (5) 
x = distance measured from the heated surface (Fig. 1) 
a = thermal diffusivity at x 0, a = (k/pc)z 
p = Density 


— INTEGRAL METHOD has been successfully developed and 
employed by Goodman! ? to solve one-dimensional heat 
conduction problems for certain boundary conditions. In the 
analysis, he approximated the temperature profile by a poly- 
nomial. The purpose of this note is to show that the same 
method can also be applied to solve the heat conduction in a slab 
of finite thickness with arbitrary heating rate and variable 
properties. Furthermore, it is also demonstrated that the 
temperature profile is better approximated by 
function than a polynomial 

The physical model is shown in Fig. 1. A slab with finite 


an exponential 


thickness, L, is initially at zero temperature. At time ¢ > 0, 
the slab is subjected to aerodynamic heating on one side and is 
insulated on the other side. The temperature u(x, f) is to be cal- 


culated. The differential equation and boundary conditions are 


as follows: 


pc(du/dt) = (0/dx)[k(Ou/dx)], Odc2z<L (1) 
u(x, 0) = 0 (2) 

k{ou(O, t)/Ox = —q (3) 

[Ou(L, t)]/Oox = 0 (4 


Using the Goodman transformation,! the above equations can be 
transformed to the corresponding simpler ones. The new vari- 


able and transformed equations are as follows: 


u = 
v= ‘ pcdu (9 
) 


0v/dt = (0/d0x)[a(dv/dx)] (6) 
wx, O) = 0 (7) 

[d2(0, t)|/Ox = —(g/a) (8) 
[dv( L, t)]/ox = 0 (9 


Integrating Eq. (6) with respect to x and using Eqs. (8) and (9), 


d f° L 
di J, vdx = q (10) 


We shall approximate the temperature profile by an exponential 


one gets 


function of the following form: 


—e gq 
vx, t) = 00,1) — L+ 
2 a 
L ( 
d (e—G/E — e~')? (11) 
21—e")a 


Substituting Eq. (11) into Eq. (10), one gets 


(d/dt)} Lv(0, t) — (L2/4)[(1 — 3e72)/(1 — e=)](g/a@)j = q (12) 


1. 10r — Exact Solution (Ref3 
1.05P >. — —Exponential Approx. eqid 
1,0. Polynomial Approx (eq tla 
= ; Fic. 2. Comparison of tem- 
%r : perature profiles. 
ar le 
a , r 
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In order to solve Eq. (12), it is necessary to know the initial 
condition 2(0, to) and tf. The fo is the time 1equired for the heat 
to penetrate from x = O0Otox = L. Both f and v(0, to) can be 
found by solving Eq. (10) with the following temperature profile 


v = 00, t)(1 — e7!)-%e7 4/2) — e-1)2 (PS) C13) 


The combination of Eqs. (8) and (13) yields 


20, to) = [(1 — e7!)/2](g/a)L (14) 


On the substitution of Eq. (13) into Eq. (10), one gets 
ly 
{ gdt = {(1 — 4e~! + 5e-2)/[4(1 — e)] f(g/a)L? (15) 
) 


Hence, one can find fp and 2(0, fo) from Eqs. (15) and (14), respec- 
tively. As soon as this is done, one can use fy and 7(0, fo) as an 
initial condition to solve Eq. (12). Finally, the temperature 
u(x, t) can be computed from Eqs. (11) and (5), consecutively. 
For constant properties and constant heat flux, the foregoing 


equations can be simplified to the following ones: 


to = {(1 — 4e7! + 5e-?)/[4(1 — e7)]f L?/a (16) 
u(O, to) = [(1 — e~')/2](q/R)L (17) 
u(O, t) = u(O, to) + (g/L)a/k)(t — to) (18) 
u(x, t) = u(0, t) — 
1—e'¢ L ( : 
tf t (e~ @/L) — @-1)2 (19) 


2 k-  2Al—e\k 
The above results have been obtained by use of an exponential 
temperature profile. If we approximate the profile by a poly- 
nomial of degree two, we shall get a set of equations somewhat 
different from the above. With a polynomial approximation, 
the resulting equations corresponding to Eqs. (16) to (19) can 


be written as follows 


to = (1/6)(L2/a) (16a) 
u(O, to) = (1/2)(g/R)L (17a) 
u(O, t) = (0, to.) + (a/R) g/L)\(t — to) (18a) 


u(x, t) = u(O, t) — (qL/k)(x/L){1 — 3/2 x/L)] (19a) 


To test the accuracy of the integral method, a numerical ex- 
ample has been carried out for the following physical values: 
L = 1/12 ft; ¢ = 1/60 hr; a = 1 ft?/hr; k = 50 Btu/hr-ft-°F; 
q = 3.6 X 10° Btu/hr-ft?. The resulting temperature profile 
is shown in Fig. 2. The exact temperature profile as calculated 
from Ref. 3 is also shown in the figure for comparison. It is 
seen from Fig. 2 that the maximum errors in temperature intro- 
duced by the polynomial and exponential approximations are 
about 6.5 percent and 3.5 percent, respectively. Hence, it can 
be concluded that the integral method with an exponential 
temperature profile can be used for heat-conduction calculation. 
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On the Entropy Change Through a Shock Layer 


James Serrin and Y. C. Whang 
Department of Mathematics and Department of Aeronautical 
Engineering, University of Minnesota, Minneapolis, Minn. 


May 11, 1961 


ws A viscous, heat-conducting gas traverses a shock 
wave its temperature, density, pressure, and entropy are 
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Fic. 1 Shock transition 
diagram. 











all increased. It might be supposed, therefore, that during the 
process of shock transition all of these quantities increase mono 
tonically from their initial to their final values. This is indeed 
true for the temperature, density, and pressure, as we shall show 
here, but it is not the case for the entropy. Rather surprisingly, 
the latter quantity first increases to a local maximum in the in- 
terior of the shock, and then undergoes a decrease to its final 
value (Fig. 2). These conclusions follow from general thermo- 
dynamic assumptions, and therefore apply to all gases of practical 
interest. 

Consider, then, a steady one-dimensional shock transition in 
an arbitrary viscous, heat-conducting fluid. As is well known, 
the transition is governed by the equations 


pu =m (m > 0) 
K(dT/dx) = m[E — '/(u — a)? — 5] (1) 
(A + 2u)(du/dx) = p + m(u — a) 


representing conservation, of mass, energy, and momentum, 

respectively. Here a, b, m are constants, \, uw, Kk are constitutive 

coefficients of the fluid (positive, but not necessarily constant), 

and the remaining notation is standard.! We require the follow- 

ing thermodynamic assumptions concerning the fluid: 

The specific heats are positive, and the following inequalities hold 
(Op/Op)s > 0, (Op/O0S)p > 0, (0?p/d7?)s > O 

where r = 1/p is the specific volume, and S is the entropy. 

The first of these inequalities simply states that the speed of 
sound is real. The second and third are slightly more special 
assumptions, but are nevertheless satisfied in all practical situa- 
tions. Now £ and p in Eqs. (1) can be expressed as functions of 
thermodynamic variables pand T. Thus Eq. (1) can be reduced 
to the following pair of equations 

K(dT/dx) = mL(u,T), (X + 2y)(du/dx) = mM(u,t) (2) 
The curves L = 0 and M = O may be plotted in the (u,T) 
plane, and have the general form shown in Fig. 1. A shock 
layer or shock transition is then a solution of Eq. (2): 
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Fic. 2. Entropy variation through shock layer (ideal gas with 
y = 7/5, (A + 2p)cp/k = 1). 
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u=ux), T=T(x), -e7<x<+« 

which, considered as a parametric curve in the (u,7) plane, joins 
the singular points Z,; and Zs. The qualitative form of the 
curves L = Oand M = 0, and the existence of a shock transition 
as shown in Fig. 1, follow from the thermodynamic assumptions 
noted above; cf. Gilbarg.? It should be observed that the shock 
curve rises to the left as x increases, because L > 0 and M < 0 
throughout the region R between the curves L = Oand V = 0. 

Consider now the variation of the thermodynamic state as one 
follows the shock curve from Z; to Z». Evidently 


dp = —mu~*du >6, dT>0O 
so that both p and T increase monotonically from Z,; to Z2. 
Moreover 


dp = (Op/Op)rdp + (O0p/0T) dT 
= (1/yT)|T(Op/dp)sdp + c,(0p/O0S)pdT| > 0 


according to our thermodynamic assumptions. It remains only 
to investigate the change in entropy. But we have 


TdS = dE + pdr = dL + Mdu 


Thus at Z; we have dS = dL/T > 0 in the direction of the shock 


>. + 


Blunt-Body Heat Transfer 


FORUM 99] 


curve, and at Z, we have dS = dL/T < 0 along the shock curve. 
In other words, although S increases as the shock begins, it must 
be decreasing over the final portion of the transition. This com- 
pletes the proof of our assertions 

For an ideal gas with constant specific heats and with Prandtl 
number (A + 2y)cp/x = 1, Eq. (2) has the solution 


Col + 3/ou? = 1/ou,2}1 + [2/(y — 1)) 173} 
as noted by Becker. In this case the entropy variation in the 
shock transition may be explicitly expressed in the form 
exp [(S — Si)/eo] = (T/T1)(p1/p)*" -= (T/T,)(u/u)? 
= tl(y — 1)/2) Mir{1 — (u?/m*)] + 1f X 


(u/u,)%~! 


5. It may be observed that 


which is plotted in Fig. 2 for y = 7/5 
in the limit !@, — 1 the maximum entropy occurs exactly when 
the velocity is half way between its initial and final values 
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sonic speed is presented. It is shown that for practical 
body shapes the viscous terms in the Navier-Stokes- 
Fourier equations are not important in the region near 
the shock, and the displacement thickness can be 
neglected. Then the flow can be approximately repre- 
sented by an inviscid flow solution having the body 
shape as boundary conditions. This solution is a good 
approximation, and is not affected by the Reynolds 
number. The flow near the wall can be analyzed by a 
boundary-layer type of analysis using appropriate 
boundary conditions. This approach permits the 
determination of the heat transfer in the region of the 
nose also at very low Reynolds numbers without the 
necessity of solving the Navier-Stokes-Fourier type of 
equations. Experimental results agree with the values 
given by the analysis. 
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SK, Honors 


and 
Awards 


of the IAS 


The IAS Honors and Awards program has for 
many years given recognition for distinguished 
achievement in the aerospace sciences. The IAS 
has a special brochure which covers the awards 
and their past recipients in detail. The following 
excerpts indicate the nature of this I[AS program. 


HONORARY FELLOWSHIP, the highest IAS honor may 
be accorded annually to two persons of eminence 
in aerospace technology, one of whom must reside 
in the U.S.A. 


IAS FELLOWS have attained a position of distinction 
in aeronautics and have made notable and valuable 
contributions in one of the aerospace sciences or 
in aerospace engineering. Not more than ten resi- 
dents of the U.S., and ten residing outside the 
U.S.A., may be elected in any one year. 


THE OCTAVE CHANUTE AWARD established by the 
IAS in 1939 honors the memory of the pioneer 
U.S. aeronautical investigator for whom it is 
named. It is given for ‘‘a notable contribution made 
by a pilot to the aerospace sciences.”’ 


THE JOHN JEFFRIES AWARD, recognizing the im- 
portance to aeronautics of scientific endeavor in 
the field of medicine, honors the memory of the 
American physician who made the earliest re- 
corded scientific observations from the air. The 
award is given for ‘outstanding contributions to 
the advancement of aeronautics through medical 
research.” 


THE LOUIS W. HILL SPACE TRANSPORTATION 
AWARD carries the largest honorarium available 
through a scientific society — $5,000 to an in- 
dividual, or up to $10,000 in the event of a ‘‘team”’ 
contribution. 


THE ROBERT M. LOSEY AWARD established in 1940 
honors the memory of Captain Robert M. Losey, a 
meteorologist, the first officer in the service of 
the United States to die in World War II. It is pre- 
sented ‘‘in recognition of outstanding contribu- 
tions tc the science of meteorology as applied to 
aeronautics.’”’ 


THE SYLVANUS ALBERT REED AWARD, endowed by 
Dr. S. A. Reed, a founder member of the IAS, to 
be presented ‘‘for a notable contribution to the 
aerospace sciences resulting from experimental or 
theoretical investigations, the beneficial influence 
of which on the development of practical aero- 
nautics is apparent.”’ 


THE LAWRENCE SPERRY AWARD, established in 
1936 by the sister and brothers of Lawrence 
Sperry, pioneer aviator and inventor who died in 
1923, is given ‘‘for a notable contribution made 
by a young man to the advancement of aero- 
nautics.”’ 


THE IAS FLIGHT TEST ENGINEERING FELLOWSHIP 
provides for two years of graduate study in flight 
test engineering leading to a master’s degree or 
doctorate in engineering. Program of study can 
be varied to meet the background and needs of 
the individual student. 





THE MINTA MARTIN NATIONAL STUDENT AWARDS, 
endowed by Glenn L. Martin, are given for out- 
standing student papers presented at the annual 
regional conferences. The regional first-prize 
papers in each division are considered for national 
awards. 
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